Nonlinear Tolerancing: Variation Simulation and Assembly Analysis with Regard to Contact Interaction of Parts

Abstract

The variation analysis is a key tool for ensuring the high quality assembly in the process of developing the technology for manufacturing of aircraft parts. One of the main factors in variations is the deviations in the positioning procedure. This paper is devoted to the development of an approach that allows taking into account the variations during positioning and merging it with the special algorithm of contact problem solving. The impact of varied boundary conditions is incorporated into an additional vector of forces that can be interpreted as reactions to the shift of supports. The obtained results are illustrated with a case of wing-to-fuselage assembly.

1. Introduction

As it was mentioned in the review by Ceglarek et al. [1], in the aircraft manufacturing and automotive industries, approximately 70% of all technological changes are carried out at the production stage, since at the design and implementation stages there is not yet enough information to accurately predict the results of technological processes. The development of digital technologies makes it possible to predict the quality of production and the probability of production defects using mathematical modeling of the assembly process. The main difficulty in predicting the outcome of a production process is the inevitable deviation of all real components in practice from the ideal constructions. Hu et al. [2] admit that the source of these assembly deviations is both manufacturing defects of individual parts and possible inaccuracies of each step of the assembly process. The cumulative effect of these deviations can lead to the fact that the assembled structure will not meet the technical requirements for reliability and functionality, as was noticed by Söderberg et al. [3]. Therefore, in the mathematical modeling of the assembly process, it is important to take into account these deviations and analyze their influence on the final deviations in the assembled structures.

The most common type of part assembly in aircraft manufacturing is riveting, which ensures the required reliability and durability of the structures. One of the important stages of riveted assembly is assembly with the use of temporary fasteners. These fasteners are installed in a specified set of assembly holes and clamp the parts together during further assembly operations.

In mass production conditions, the arrangement of temporary fasteners cannot be adjusted to the individual characteristics of incoming parts. Therefore, a single template for installing temporary fasteners is developed for each joint. Despite possible random assembly deviations, this template should eliminate the gap between parts. As admitted by Weber [4], the development of effective templates for temporary fasteners is one of the actively used methods for optimizing the assembly process in aircraft manufacturing.

The main approach to handling problems of analyzing the assembly process, taking into account the deviations, is statistical modeling. The general solution scheme in this case is based on the Monte Carlo method and consists of cycled modeling of the outcome of the assembly process for random initial deviations. Based on the data obtained, the required probabilistic characteristics of the analyzed process are determined using statistical methods. It is necessary to have a large set of implementations of assembly deviations to carry out the statistical analysis. Additionally, since most of the parts used in aircraft manufacturing are compliant structures, it is important to take into account the deformation of the parts during their contact interaction.

The standard method that allows taking into account assembly deviations is the direct Monte Carlo method, as described by Gao et al. in [5]. The realizations of the initial deviations are obtained by generating random numbers from given distributions, and the modeling of the assembly process is carried out using the classical finite element method. However, the use of the finite element method in the series of calculations for large and complex structures is a very labor-intensive approach.

In 1997, Liu and Hu [6] proposed the method of influence coefficients (MIC), which makes it possible to establish a linear relationship between the initial deviations and the characteristics of the analyzed assembly process. This method significantly reduces the time required for statistical modeling. In this case, however, it is assumed that contact between the assembled structural elements occurs only at predetermined points (for example, at welding points), which makes MIC a purely linear method of analysis and significantly reduces the reliability of results in many application areas. However, the simplicity and speed of calculations have made MIC an extremely popular research method in many fields, including aircraft manufacturing [7,8,9,10,11,12].

In 2008, Warmefjord et al. [13] proposed the development of MIC by taking into account the possibility of contact between assembled structural elements at arbitrary points in the junction area. In this work, the method of influence coefficients is combined with methods of direct search for contact points. This approach is quite widely used. For example, articles [3,14,15,16,17,18] describe MIC application to the analysis of various assembly processes. However, this approach is difficult to implement and does not allow working with sufficiently detailed models.

In the approach proposed in 2010 by Lupuleac et al. in [19], the variational formulation of the contact problem (e.g., see the book of Wriggers [20]) is reduced to a quadratic programming problem (QPP). The most common methods for QPP were investigated, and specialized solvers have been developed for handling problem features. The dual active set method proposed by Goldfarb [21] is proven to be fast and efficient for small to medium-sized problems (up to ten thousand unknowns, see Burton and Toint [22]). The problem requires focusing on the modification made by Powell [23] that is preferable for ill-conditioned problems. The warm start is easy to implement for active set methods, and one can make use of previously found solutions for subsequent computations. The next one, the interior-point method, is a polynomial-time algorithm that allows for solving a wide range of optimization problems, such as linear, convex quadratic programming, etc. [24]. This method is known to be one of the main methods used for solving large sparse problems and is implemented in many commercial optimization software packages (CPLEX, IMSL, MATLAB, and others). The appropriate choice of a starting point and the advances in the solution of the linear system for determining the Newton direction make the interior-point method applicable to the types of problems considered. Reformulating the QPP in a dual and special “relative” form [25] provides the simple box constraints and leads to the Newton projection method [26], which has a quadratic convergence rate. As a result, the contact problem for deformable parts is solved quickly and without loss of accuracy. This methodology was significantly expanded in subsequent works devoted to modeling the assembly process in the aircraft industry [27,28,29,30].

This study attempts to combine the approaches proposed in [9] and [19]. Thus, the approach presented here combines the simplicity and speed of MIC with the accuracy of contact analysis.

Note that the methods taking into account the deviations of parts are widely used for the optimization and development of assembly technology for aircraft structures. In particular, the optimization of the wing-to-fuselage assembly process for Airbus A350 aircraft is considered in [27]. The optimization of the wing assembly process for Airbus A320 is discussed in [28]. The procedure for optimizing the S19 splice joint assembly in the tail section of A350 is described in [29].

The remainder of the article is organized as follows: Section 2 is devoted to the mathematical approaches used in the work. Section 2.1 and Section 2.2 describe a methodology that allows positioning deviations to be easily taken into account when solving contact problems in assembly modeling for compliant structures. Section 2.3 describes special methods that allow one to quickly solve the numerous quadratic programming problems that arise when modeling assembly processes. Section 2.4 discusses the statistical analysis of the results. Section 3 describes the results of mathematical modeling. The surrogate model of the wing-to-fuselage upper joint that is used for the numerical experiments is presented in Section 3.1. Section 3.2 illustrates the influence of deviations in the positions of support points on the gap between assembled parts. The choice of the most appropriate numerical method for the considered problem is discussed in Section 3.3. The main modeling results on the influence of positioning deviations on the assembly quality and resulting deviations are provided in Section 3.4. The emulation of positioning deviations by setting the initial gaps between assembled parts is discussed in Section 3.5. The main conclusions of the work are contained in Section 4.

2. Methodology Concept

2.1. Problem Statement

This section is devoted to the problem statement and description of the numerical procedure.

Let us consider the finite element model (FEM) of two panels to be assembled, as shown in Figure 1a. The parts may come into contact in the area where finite element nodes are marked by colored balls; further on, this area is referred to as the junction area. The panels are fixed along the edges (clamps are drawn in Figure 1a). According to the FEM, it is possible to split the nodes into three groups: nodes that can be brought into contact—the computational ones; nodes where fixations are set—nodes to be varied and all the rest. The degrees of freedom (DOF) in computational nodes are referred to as $U_C={\left(U_{Ci}\right)}_{i=1}^{N_C}$, the DOF in nodes to be varied—as $U_V={\left(U_{Vi}\right)}_{i=1}^{N_V}$, the DOF in the rest of nodes—as $U_R={\left(U_{Ri}\right)}_{i=1}^{N_R}$. All the groups of DOF are represented in Figure 1b.

Thus, the vector of unknown displacements in FE nodes $U={\left(U_i\right)}_{i=1}^N$ can be constructed:

$$U={\left(U_C,U_V,U_R\right)}^T,\hspace{1em}N=N_C+N_V+N_R$$

(1)

If the possibility of contact interaction is not taken into account, then the vector $U$ is the solution of the linear system:

$$K\cdot U=F,$$

(2)

where $K={\left(k_{ij}\right)}_{i=1,j=1}^{N,N}$ is the global stiffness matrix of the FEM, $F={\left(f_i\right)}_{i=1}^N$ is the vector of applied loads that includes both forces and moments. One can note that system (2) is obtained by discretization of the equations of elasticity theory. Since the assembled parts are mainly flexible panels, shell theory is usually used to describe the stress-strain state of the parts.

The boundary conditions correspond to the fixations of parts on the assembly stand, for example, with the help of clamps, as shown in Figure 1a:

$$U_V=U_{V0},$$

(3)

where $U_{V0}\in R^l$ is the vector of fixed degrees of freedoms.

The contact interaction between the considered parts is modeled using the node-to-node contact model (see Figure 1b). According to this model, only predefined pairs of nodes can come into contact. Note that the use of such a model for this kind of problem is justified since the installed fasteners prevent significant tangential displacements of parts relative to each other.

Using the node-to-node contact model, the nonpenetration conditions can be formulated as:

$$A\cdot U\le G,$$

(4)

where $A\in R^{M\times N}$ is a linear operator that defines nonpenetration between the pairs of nodes in the junction area, $G\in R^M$ is the vector of initial gaps between the panels. It can be mentioned that $A=0$ outside the junction area, so it is possible to rewrite (4):

$$A_C\cdot U_C\le G,$$

(5)

Summing up the relations (2)–(5), the problem of calculating vector $U$ can be stated in a variational form as the minimization problem:

$$\begin{array}{c}\underset{U\in S}{\min }\left(\frac{1}{2}{U}^T\cdot K\cdot U-F^T\cdot U\right)\\ S=\left\{\left.U\right|A_C\cdot U_C\le G,U_V=U_{V0}\right\},\end{array}$$

(6)

The statement (6) of the contact problem includes all the degrees of freedom of the structure taken into account in the FEM, which dimension $N$ can reach ${10}^6$–${10}^7$ variables for industrial applications. This factor makes the solution of the contact problem in statement (6) quite resource-intensive, especially if multiple calculations are required.

On the contrary, the number of DOF in the junction area $N_C$ is much smaller than the dimension $N$. The static condensation technique [31] can be applied in order to reduce the problem dimension. This approach is described in the next subsection.

2.2. Problem Reformulation

Let us divide the problem (6) into blocks according to the structure of the displacement vector $U$ from (1):

$$K=\left(\begin{array}{c}\begin{array}{cc}{K}_{CC}& K_{CV}{K}_{CR}\\ K_{CV}^T& K_{VV}{K}_{VR}\end{array}\\ \begin{array}{ccc}{K}_{CR}^T& K_{VR}^T& K_{RR}\end{array}\end{array}\right),F=\left(\begin{array}{c}\begin{array}{c}{F}_C\\ F_V\end{array}\\ F_R\end{array}\right)$$

(7)

The loads are applied only inside the junction area, so $F_V=F_R=0$. The system (2) providing a global minimum of (6) without accounting for $S$ is rewritten:

$$\left(\begin{array}{c}\begin{array}{cc}{K}_{CC}& K_{CV}{K}_{CR}\\ K_{CV}^T& K_{VV}{K}_{VR}\end{array}\\ \begin{array}{ccc}{K}_{CR}^T& K_{VR}^T& K_{RR}\end{array}\end{array}\right)\cdot \left(\begin{array}{c}\begin{array}{c}{U}_C\\ U_V\end{array}\\ U_R\end{array}\right)=\left(\begin{array}{c}\begin{array}{c}{F}_C\\ 0\end{array}\\ 0\end{array}\right)$$

(8)

Opening the brackets, (8) transforms into:

$$\begin{array}{c}{K}_{CC}\cdot U_C+K_{CV}\cdot U_V+K_{CR}\cdot U_R=F_C,\\ K_{CV}^T\cdot U_C+K_{VV}\cdot U_V+K_{VR}\cdot U_R=0,\\ K_{CR}^T\cdot U_C+K_{VR}^T\cdot U_V+K_{RR}\cdot U_R=0\end{array}$$

(9)

Third relation provides the formula for $U_R$:

$$U_R=-K_{RR}^{-1}\cdot \left(K_{CR}^T\cdot U_C+K_{VR}^T\cdot U_V\right)$$

(10)

Then, $U_R$ from (10) is substituted into the first relation of system (9):

$$K_{CC}\cdot U_C-K_{CR}\cdot K_{RR}^{-1}\cdot \left(K_{CR}^T\cdot U_C+K_{VR}^T\cdot U_V\right)=F_C-K_{CV}\cdot U_V$$

(11)

It is recalled that $U_V=U_{V0}$ and (11) is simplified:

$$\left(K_{CC}-K_{CR}\cdot K_{RR}^{-1}\cdot K_{CR}^T\right)\cdot U_C=F_C-K_{CV}\cdot U_{V0}+K_{CR}\cdot K_{RR}^{-1}\cdot K_{VR}^T\cdot U_{V0}$$

(12)

From (12), it is possible to derive matrix $K_C=K_{CC}-K_{CR}\cdot K_{RR}^{-1}\cdot K_{CR}^T$ that can be interpreted as the reduced stiffness matrix and the force vector for variations of boundary conditions:

$$F_{var}=\left(-K_{CV}+K_{CR}\cdot K_{RR}^{-1}\cdot K_{VR}^T\right)\cdot U_{V0}=M_{var}\cdot U_{V0}$$

(13)

Adapting the notations from [6], $M_{var}=\left(-K_{CV}+K_{CR}\cdot K_{RR}^{-1}\cdot K_{VR}^T\right)$ can be called the matrix of influence coefficients.

Given all the previous transformations the reduced minimization problem equivalent to (6) can be formulated:

$$\begin{array}{c}\underset{U_C\in S_C}{\min }\left(\frac{1}{2}{U}_C^T\cdot K_C\cdot U_C-\left(F_C^T+{\left(M_{var}\cdot U_{V0}\right)}^T\right)\cdot U_C\right)\\ S_C=\left\{\left.U_C\right|A_C\cdot U_C\le G\right\},\end{array}$$

(14)

One can see that matrix $K_C$ can be computed once, the same is applicable for matrix $M_{var}$. A series of assembly simulation problems can then be constructed by only using the input vector of loads $F_C$ in the junction area, the vector of initial gaps $G$ and the vector of boundary values $U_{V0}$.

In the case of matrix $K$, its blocks from (7) are not known; the static finite element analysis (FEA) can be exploited for calculating matrix $K_C$ and $M_{var}$.

Let us suppose that $U_{V0}=0$, and ${\left(F_C\right)}_j=\left\{\begin{array}{c}1,\hspace{1em}j=i\\ 0,\hspace{1em}j\ne i\end{array},j\in \left[1,N_C\right]\right.$ for i-th DOF in the junction area. The FEA analysis is performed, and the corresponding $U_C$ is computed. Then, it can be concluded from (12) that:

$$U_C={K_C^{-1}\cdot F}_C={\left\{K_{Cj}^{-1}\right\}}_{j=1}^{N_C}$$

(15)

After collecting all the columns of $K_C^{-1}$ the matrix inversion is required to obtain $K_C$ itself.

The columns of the matrix $M_{var}$ are obtained by setting ${\left(U_{V_0}\right)}_j=\left\{\begin{array}{c}1,\hspace{1em}j=i\\ 0,\hspace{1em}j\ne i\end{array}\right.$ for the i-th DOF in the nodes to be varied ($j\in \left[1,N_V\right]$) and $F_C=0$. $U_C$ is computed by static FEA and revealing (12):

$$K_C\cdot U_C={\left\{M_{varj}\right\}}_{j=1}^{N_C},$$

(16)

So, the i-th column of $M_{var}$ equals $K_C\cdot U_C$.

2.3. Algorithms for Quadratic Programming Problems

Obtaining matrices $K_C$ and $M_{var}$ for the problem (14), the variational analysis comes down to solving a set of quadratic programming problems with different input data. To achieve significant performance improvements over standard quadratic programming tools, the specific features that characterize the considered problem can be used together with:

• reformulating the quadratic programming problem in order to reduce the number of unknowns and simplify the constraints;
• using specialized optimization algorithms adapted to assembly problems;
• dividing the calculation process into a preprocessing stage (time-consuming one) and a gap calculation stage (fast one).

2.3.1. Reformulation of the Quadratic Programming Problem

In assembly problems, the formulation of the original quadratic programming problem in the form most preferable for the optimization method can help reduce computation time. More specifically, this formulation should have a low dimension with a simple constraint structure. Let us consider two additional formulations of the minimization problem (14) that meet the above requirements. The dual formulation of problem (14) is given by:

$$\underset{{\mathsf{\Lambda }}_C>0}{\max }\left(-\frac{1}{2}{\mathsf{\Lambda }}_C^T\cdot Q_C\cdot {\mathsf{\Lambda }}_C+P_C^T\cdot {\mathsf{\Lambda }}_C+H\right),$$

(17)

where ${\mathsf{\Lambda }}_C\in R^M$ is the vector of Lagrange multipliers, $Q_C={A_C}^T{K}_C^{-1}{A}_C\in R^{M\times M}$ is a symmetric positive definite dense matrix, $P_C=A_C^T\cdot K_C^{-1}\cdot T_C-G\in R^M$, $H=-\frac{1}{2}{T}_C\cdot K_C^{-1}{\cdot T}_C\in R^1$, $T_C=F_C^{}+M_{var}\cdot U_{V0}$. The physical meaning of the unknown vector ${\mathsf{\Lambda }}_C$ is the reaction forces arising in the contact nodes.

The relative formulation [25] of problem (14) is written as:

$$\begin{array}{c}\underset{D_C\in {\tilde{S}}_C}{\min }\left(\frac{1}{2}{D}_C^T\cdot {\tilde{K}}_C\cdot D_C-{\tilde{F}}_C\cdot D_C\right)\\ {\tilde{S}}_C=\left\{\left.D_C\right|D_C\le G\right\},\end{array}$$

(18)

where $D_C=A_C\cdot U_C\in R^M$ is a vector of relative displacements, ${\tilde{K}}_C=Q_C^{-1}\in R^{M\times M}$ is a symmetric positive definite dense matrix, ${\tilde{F}}_C={\tilde{K}}_C{\cdot A}_C^T{\cdot K}_C^{-1}{T}_C\in R^M$.

Solving problem (17) and (18) is equivalent to solving (14). The solution of (14) $U_C$ is expressed from the solutions of (17) ${\mathsf{\Lambda }}_C$ and (18) $D_C$ using Formula (19):

$$\begin{array}{c}{U}_C=K_C^{-1}\left(T_C-{A_C\cdot \Lambda }_C\right),\\ {\mathsf{\Lambda }}_C={\tilde{F}}_C-{\tilde{K}}_C{\cdot D}_C\end{array}$$

(19)

The formulations (17) and (18) reduce the dimension of the quadratic programming problem (14) from $N$ unknowns to $M$ ($N\ge M$), the number of constraints remains unchanged, but the linear constraints $S_C$ in general form are replaced by the bound constraints ${\mathsf{\Lambda }}_C>0$ and ${\tilde{S}}_C$, respectively.

2.3.2. Adaptation of Algorithms

Problem (14) itself is a convex quadratic programming problem, the solution of which exists and is unique. The optimization methods are faced with the challenge of not only finding a solution with some predetermined accuracy but also finding it quickly since it is necessary to solve a large number of similar problems. The dual Goldfarb–Idnani active set method [21,23], a Newton projection method [26], the primal-dual interior-point method [32], and Lemke’s method [33] from this point of view are discussed further.

Depending on the specifics of the optimization problem, the methods are distinct because of different basic ideas. Features of assembly problems lead to the formulation of a quadratic programming problem with a block-diagonal structure of the matrix $K_C$, a sparse structure of the constraint matrix $A_C$, ill-conditioned, symmetric, and positive-definite Hessian matrices $K_C$, $Q_C,\mathrm{and}$ ${\tilde{K}}_C$. In addition, variation simulation provides a set of quadratic programming problems where the Hessian matrix and constraint matrix are the same, but the force and gap vectors are different. These features were used to adapt the listed methods to effectively solve assembly problems. The Newton projection method and the primal-dual interior-point method have an iterative structure, where at each iteration, the system of linear equations is solved. The adaptation (see

Table 1. Adaptation of optimization methods to assembly problems.

Newton Projection Method	Dual Goldfarb–Idnani Active Set Method	Primal-Dual Interior-Point Method	Lemke’s Method
A method for recalculating constraints to reduce the number of iterations of the method. A method for solving a system of equations based on a combination of direct and iterative approaches. A modification of direct methods for solving the system based on the Sherman–Morrison formula. A method for finding step length based on the golden section method. Using a warm start.	Taking into account the sparse structure of the constraint matrix $A_C$. Taking into account the block-diagonal structure of the Hessian matrix $K_C$. Improving the procedure for adding a new constraint to the working set, based on the physical meaning of the problem. Using a warm start.	A method for searching for a “feasible” starting point, taking into account the specifics of assembly problems. The method makes it possible to reduce the number of iterations. Proposing a preconditioner for solving a system of linear equations reducing iteration process of the system solution. A modification of direct methods for solving the system based on the Sherman–Morrison formula.	Taking into account the sparse structure of the constraint matrix $A_C$. Taking into account the block-diagonal structure of the Hessian matrix $K_C$.

) of these methods includes improving system solving (choice between direct and iterative methods, selection of preconditioner for iterative approaches) and reducing the number of iterations (choice of starting point, strategy for choosing a new approximation at each iteration). The dual Goldfarb–Idnani active set method and Lemke’s method are based on looking over the set of active constraints. Adaptation was based on finding a strategy for modifying the working set of constraints. For all optimization methods, the structure of the Hessian and constraint matrix was taken into account to improve memory management and reduce the number of actions for matrix-vector operations. All methods presented in Table 1 are implemented in the computer code and included in the ASRP software package [30].

Note that the optimization problem under consideration is strictly convex (for all presented formulations) and therefore has a unique solution, which is found by all the algorithms presented here. The exit condition for all algorithms is the fulfillment of optimality conditions with an accuracy of 10⁻⁷.

2.3.3. Preprocessing Stage and a Gap Calculation Stage

Preparing matrices $Q_C$, ${\tilde{K}}_C$, and vectors $P_C$, ${\tilde{F}}_C$ for problems (17) and (18) requires additional computational costs. However, since the matrices $K_C$ and $A_C$ are the same for a set of problems formed by variation analysis, the calculation of the necessary auxiliary matrices can be carried out once for each assembly model at the preprocessing stage of model preparation. Taking into account the sparse structure of the matrix $A_C$ and the block-diagonal structure of the matrix $K_C$ allows one to reduce the computation time of the preprocessing stage. Preparing the relative formulation (18) is the most labor-intensive. The ill-conditioning of the stiffness matrix $K_C$ does not allow one to obtain the matrices ${\tilde{K}}_C={\left({A_C}^T\cdot K_C^{-1}\cdot A_C\right)}^{-1}$ and ${\tilde{K}}_C\cdot A_C^T{\cdot K}_C^{-1}$ with sufficient accuracy, when using standard inversion methods, due to the accumulation of rounding errors, and forces one to use more accurate methods of matrix inversion and multiplication [25].

Another part of data preparation is directly related to the specifics of optimization methods. Some operations, such as Hessian inversion, replacement of variables, construction of auxiliary matrices for the formation of systems of equations, and associated preconditioners necessary for the operation of the methods, can be performed in advance. The set of necessary auxiliary operations depends on the method used.

By moving all labor-intensive operations beyond the immediate solution of the quadratic programming problem, it allows one to significantly speed up the solution process at the gap calculation stage.

2.4. Methodology Application

Given the possibility of making fast multiple computations with the algorithms described in the previous section, the methodology of assembly analysis can be derived that accounts for variations from different sources.

Firstly, varying the force vector $F_C$ means that different fastening configurations are considered. Different vectors $G$ represent the influence of the initial gap on the assembly [34]. Varying $U_{V_0}$ corresponds to changes in the initial positioning of the parts before the assembly.

All the variations are transferred to the preprocessing stage of the QP problem, where the most matrix computations are performed (see Figure 2). Then, the set of corresponding entities is constructed and solved. Then, the results are collected and analyzed.

In order to compare simulation results, the quality criterion is to be introduced. Let us denote $U_{res}$ the solution of problem (8), $\delta$ is the threshold gap value. If the resulting gap $G_{res}=G-A_c\cdot U_{res}$ undergoes $\delta$, it is assumed that the contact is achieved. The defect gap vector $d(\delta )\in R^k$ can be defined as follows for each DOF in the junction area:

$$d_j(\delta )=\left\{\begin{array}{ll}1,& {\left(G_{res}\right)}_j\ge \delta \\ 0,& {\left(G_{res}\right)}_j<\delta \end{array}\right.$$

(20)

The averaged sum of the $d(\delta )$ components shows how the gap is eliminated for the given data case $G,F_C,{\mathrm{and}U}_{V_0}$. Further on, it is called the probability of defect:

$$P(\delta )=\frac{1}{k}\sum _{j=1}^k{d}_j(\delta )$$

(21)

The lower the probability of defect $P(\delta )$, the better the quality of the assembly.

Using this methodology, it is possible to simulate what-if scenarios when it is necessary to examine the quality of the assembly under the conditions derived from the data provided by the assembly line (e.g., the information on initial gaps between parts measured for the previously manufactured aircrafts).

Another possibility is to check the given assembly against the unknown boundary conditions a priori. Typical data known for part positioning is tolerance $T$ at each boundary point (it represents how far the actual position can vary from its nominal location): ${\left(U_{V_0}\right)}_j\in \left[-T_j,T_j\right]$ for each boundary node to be varied $j\in \left[1,l\right]$. The statistical distribution for the variation of each point ${\left(U_{V_0}\right)}_j$ is typically assumed to be normal (Gaussian) with a zero mean ${\mu }_j=0$ and with a standard deviation ${\sigma }_j=T_j/3$ [35]. Thus, $U_{V_0}$ is constructed from randomly generated ${\left(U_{V_0}\right)}_j$ by sampling corresponding normal distributions.

It is worth investigating how the variations $U_{V_0}$ affect the resulting displacements after the contact interaction. In [6], the same research was conducted, and the linear dependence was obtained. Since the approach presented in this paper is nonlinear, the dependence between source variations and variations in the points of observation is more complex. At the same time, using the above-described computation techniques allows for a speed of calculations that is comparable to that of linear models.

3. Results

3.1. Assembly Model of Wing-to-Fuselage Upper Joint

The described approach is applied to the model of aircraft wing-to-fuselage assembly. Figure 3a shows the scheme of the upper joint of the wing (gray) to a part of the fuselage called the cruciform (green). In the top view in Figure 3b, dots show the computational nodes that form the junction area (the area of possible contact), and gray circles show the holes where fasteners can be installed. The assembly model is full-scale. The size of the junction area is approximately 6 m by 0.3 m. There are 3612 computational nodes and 256 holes in the junction area. This model does not correspond to any specific type of aircraft, but qualitatively reproduces the characteristics of the real assembly process. The junction area is colored according to the gap between the parts.

The three edges of the cruciform shown in red in Figure 3a are attached to the central wing box, which is very rigid. Therefore, in the finite element model presented in Figure 4, all displacements and rotations in the nodes located on these edges are restricted. The wing is fixed at four support points depicted in Figure 3 and Figure 4 with red triangles. At these four points, all displacements except the vertical ones are restricted, and the vertical displacements are to be varied. A finite element mesh with these boundary conditions is used to calculate the reduced stiffness matrix $K_C$ and matrix of influence coefficients $M_{var}$ of the assembly, as described in Section 2.2.

3.2. Examples of the Influence Deviations in the Positions of Support Points on the Gap between Assembled Parts

In all further numerical experiments in this section, it is assumed that the parts are perfectly adjacent to each other: initial gap between parts is equal to zero with no deviations of support positions. So, the source of the gap deviations in the junction area is concentrated at the variation of the vertical positions of the support points of the wing.

Note that these variations are nonlinearly reflected in the gap at the junction area. For example, Figure 5a,b show the results of moving all support points in vertical directions up and down, respectively. When all the support points deflect 3 mm downwards, the gap is opened for approximately 3 mm, and when moving support points 3 mm upwards, due to contact interaction, a gap opens up to approximately 0.5 mm at the edge of the wing.

It is important to evaluate the impact of random deviations in the positions of the support points on the quality of the assembly. These random deviations cause various residual gaps between parts. The initial assembly stage consists of installing approximately 10% of all fasteners. As a rule, these fasteners are installed manually, which takes a relatively long time and incurs labor costs. The technological requirement is that all residual gaps with these fasteners installed are less than the specified value. The correct arrangement of fasteners at the initial assembly stage is a very important task in aircraft assembly. In this regard, the calculations of residual gaps between parts with 10% installed fasteners are considered in this work as the benchmarks. Figure 6 shows examples of such gaps without fasteners (a) and with 26 fasteners installed (that is approximately 10% of all fasteners), which corresponds to the end of the initial assembly stage (b).

3.3. Choosing a Method for Solving a Quadratic Programming Problem to Reduce Computation Time

Statistical analysis involves a lot of calculations, so choosing the fastest method is very important. Let us consider the computation time using the example of the gaps shown in Figure 6.

Figure 7 presents the computation time of the gap calculation stage spent by different optimization methods used to solve (14), (17), and (18). Similar trends are observed for different types of gaps caused by random deviations in the position of wing support. The fastest methods are the dual Goldfarb–Idnani active set method and the Newton projection method, used for the relative formulation (18). The computation time estimate of the optimization methods depends on the number of active constraints at optimum points, which is related to the number of installed fasteners (see Figure 8a). In particular, the dual Goldfarb–Idnani active set method has a proportional and Newton projection method has an inverse proportional dependence of computation time on the number of active constraints at optimum for relative problem formulation. Figure 8b helps to choose between the methods based on the number of fasteners to be installed. That is, with a small number of fixtures (less than 10%), the dual Goldfarb–Idnani active set method is preferable to use; otherwise, the Newton projection method is faster.

Note that despite the fact that the gaps for Cases 2 and 4 (see Figure 6) look qualitatively different, the dependence of the active constraints on the number of installed fasteners for these cases is almost the same (Figure 8a). This is due to the fact that both of these gaps are positive throughout the junction area, and contact occurs only after the installation of fasteners. For Cases 1 and 3, in some areas, contact is achieved before fasteners are installed. This is illustrated in Figure 9a, where the contact nodes are colored pink. From Figure 9b, it can be seen that with 10% fasteners installed, the contact nodes for Cases 2 and 4 are located only near the fasteners, while for Cases 1 and 3, they are also present in the middle of the junction area.

Further calculations are carried out without fasteners and with 10% fasteners installed; therefore, the relative formulation of the QP problem is solved by the ASM, since for these cases it is the fastest (see Figure 7 and Figure 8b).

3.4. The Influence of Random Deviations in the Positions of Support Points on a Assembly Process

Figure 10 shows examples of statistical distributions of support point deviations for different sample sizes. The column width in the histograms is 0.5 mm. The Gaussian distribution with a zero mean ${\mu }_j=0$ and standard deviation ${\sigma }_j=3$ mm is used.

Figure 11 shows the statistical curves representing the probability of defect $P(\delta )$ for different $\delta$. The $\delta$ values in mm are plotted on the horizontal axis, and the vertical axis shows the $P(\delta )$ as a percentage, calculated from the entire sample of the support point deviations.

The statistical curves without fasteners are slightly different (see Figure 11a) for different sample sizes, but after the fastener installation, this difference ceases to be noticeable starting from a sample size of 200 (Figure 11b).

Figure 12 illustrates the influence of random deviations of support points on the distribution of normal displacements of some nodes in the junction area. No fasteners are installed, but the wing and cruciform parts experience contact interaction. Figure 12a,b correspond to the case of deviation of one support point and Figure 12c to the case when two support points deviate.

As can be seen in Figure 12, despite the fact that the source variations are given by a Gaussian distribution, the deviations at the observation points have a completely different behavior. This is due to the fact that contact interactions are taken into account, which leads to the nonlinearity of the developed tolerance model. The linear model can change the parameters, but cannot change the structure of distributions (compare to Figure 8 and Figure 9 in work [6]).

3.5. Comparison of Gaps Caused by Deviations of Support Points with Generated Initial Gaps

The vector of nonuniform initial gaps between the assembled parts $G$ can be used as an alternative source of deviations in tolerancing and variational analysis. It is of interest to compare two approaches to generating deviations during the assembly process. This comparison is the subject of the present section.

The initial gap between the assembled parts can be modeled as the homogeneous Gaussian random field $G\left(x,y\right)$ with mean ${\mu }_{RF}$, standard deviation ${\sigma }_{RF}$, and exponential correlation function $\rho (x,y;{\alpha }_{RF})=\exp \left(-{\alpha }_{RF}^2\left(x^2+y^2\right)/2\right)$, where $(x,y)$ is the local coordinate in the junction area [27]. The mean ${\mu }_{RF}$ represents the constant value of the initial gap in the junction area and is assumed to be zero. Thus, the modeled initial gap is defined by two parameters: the standard deviation ${\sigma }_{RF}$ (which corresponds to the amplitude of the initial gap) and the correlation coefficient ${\alpha }_{RF}$ (which corresponds to the curvature of the initial gap) [28].

The parameters of the modeled initial gaps ${\sigma }_{RF}$ and $\alpha$ were chosen in order to replicate the deviations in the vertical positions of support points. Figure 13 shows the result of statistical analysis for two groups consisting of 5000 samples of residual gaps each. The first group was obtained by the variation of vertical positions for all four support points (${\mu }_j=0$, ${\sigma }_j=3$ mm) and the second group was built on the base of the generated initial gaps (${\mu }_{RF}=0$, ${\sigma }_{RF}=2$ mm, ${\alpha }_{RF}=0.000145$). As can be seen from Figure 13, the statistical curves for the two groups are quite close to each other, especially in the case of 10% of installed fasteners.

The same conclusion can be drawn from Figure 14, which shows four examples of the calculated residual gaps (with and without installed fasteners) based on the generated initial gaps. Comparing it with Figure 6, showing the gaps obtained by the deviation of support positions, one can see that the resulting distributions are very similar.

Figure 15 presents the comparison distributions of maximal and mean values for the same groups of residual gap samples. As can be seen, the distributions of the mean and maximum gaps are quite close, both qualitatively and quantitatively. Therefore, the deviations in the shape of the assembled parts can be emulated with sufficient accuracy by setting the initial gap with proper settings.

Figure 16 provides the so-called local statistics, which are the probability of obtaining a residual gap at a given point in the junction area within the prescribed range. The same two groups of gap samples are used. The distribution of local statistics provides important information for the assembly engineer, since it characterizes the quality of contact between the parts being connected at a given stage of assembly with regard to deviations of parts. In addition, local statistics make it possible to find weak points in the assembly technology under study, showing in which places the arrangement of fasteners may need to be changed.

Based on Figure 15 and Figure 16, one can draw similar conclusions about the proximity of the final gaps obtained by different methods of generating deviations.

The analysis carried out in this section shows that positioning deviations can be adequately emulated by setting the initial gaps between assembled parts. However, selecting parameters for these initial gaps is a nontrivial task, the solution of which is impossible in the absence of information about the effect of positioning deviations.

4. Conclusions

The presented work continues the line of research started by the famous work of Liu and Hu [6], where the method of influence coefficients was proposed. The proposed generalization makes it possible to quickly and easily evaluate deviations in the resulting assembly and the quality of contact based on positioning deviations. At the same time, unlike classical MIC, the proposed method provides a mathematically correct solution to the contact problem that inevitably arises when modeling the assembly process.

Note that, along with positioning deviations, the quality of the final product is influenced by deviations in the assembled parts that appeared as a result of previous assembly operations. Such deviations can be modeled by setting a nonzero initial gap between parts, in particular according to the methods proposed in [34]. For building the most complete deviation model, these two approaches can be combined.

The main novelty of the proposed method is the full implementation of deviation modeling in the process of solving the contact problem. In this context, taking into account positioning deviations is not fundamentally different from taking into account the action of forces (for example, from fastening elements). The proposed technique allows for nonlinear tolerancing with a computation speed that is comparable to the speed of solving linear tolerancing problems of similar dimensions.

Note that the developed modeling methodology is not limited to the aerospace industry. It can be used in any industry where it is necessary to assemble structures from large, compliant parts. In particular, this takes place in the automotive industry, the assembly of cars for high-speed trains, shipbuilding, and other industries.