Numerical Approach for Detecting the Resonance Effects of Drilling during Assembly of Aircraft Structures †

: This paper is devoted to the development of a numerical approach that allows quick de-tection of the conditions favorable for the beginning of noticeable vibrations during drilling. The main novelty of the proposed approach lies in taking into account the deviations of the assembled compliant parts during non-stationary contact analysis by means of variation simulation. The ap-proaches to stationary analysis of assembly quality are expanded and generalized for modeling such non-stationary effects as vibration and resonance. The numerical procedure is based on modeling the stress–strain state of the assembled structures by solving the corresponding transient contact problem. The use of Guyan reduction, the node-to-node contact model and the application of the generalized α method allow the reformulation of the contact problem in terms of a series of quadratic programming problems. The algorithm is thoroughly tested and validated with commercial software. The efficiency of the developed numerical procedure is illustrated by analysis of the test joints of two aircraft panels. The unsteady process of drilling the panels with periodic drilling force was simulated. The influence of deviations in the shape of the parts on the non-stationary interlayer gap was modeled by setting different initial gaps between parts. It is shown that the oscillation amplitudes of the interlayer gap depend on the initial gaps and do not correlate with the mean value of the stationary residual gap. Thus, non-stationary analysis provides new information about the quality of the assembly process, and it should be applied if the assembly process includes periodic impact on the assembled parts.


Introduction
The assembly of a commercial aircraft involves multiple drilling and reaming operations in order to make holes for fastener installation [1,2]. A common hole-producing operation in aircraft assembly is drilling holes through multiple layers of sheet metal or composite materials in order to fasten them together [3,4]. One of the most important conditions to be met during drilling is to ensure that the interlayer gap between the assembled parts is less than the specified value (usually a few tenths of a millimeter). Violation of this condition can lead to problems such as burr formation, ingress of chips between the assembled parts and misalignment of the drilled holes [5][6][7][8].
Nowadays, mathematical modeling is increasingly used to optimize the assembly process [9][10][11]. In particular, the drill loads are taken into account in analysis and optimization of the assembly process [12]. Luo et al. [13] investigated the drilling parameters of Al7075-T6 aerospace aluminum alloy and the chip formation using finite element analysis. Liu et al. [14] optimized the drilling process to reduce the interlayer gap through numerical simulation and experimentation. Tinkloh et al. [15] investigated the effect of composite waviness on the drilling process.
There are a lot of studies devoted to vibration during drilling. Savilov et al. [16] proposed vibration suppression methods. Zhang et al. [17] developed an ambient vibration control system to improve the precision and efficiency of robotic drilling systems. In addition, methods for improving the drilling process using vibrations have been actively developed recently [18][19][20][21]. Therefore, non-stationary modeling of the assembly process, taking into account deviations of the parts, becomes very relevant. At the same time, the software tools for simulation and optimization of the assembly process use stationary contact analysis coupled with variation and statistical analysis [22], and the drilling loads are set as constant. This excludes from consideration the simulation of such phenomena as vibration and resonance, which can occur during drilling. In this paper, an attempt is made to expand the existing methodology for modeling these non-stationary phenomena.
In the numerical analysis of the assembly, it is necessary to take into account the shape variations of the assembled parts caused by manufacturing discrepancies, fixation tolerances, etc. The accounting for compliance of the parts can be performed with finite element modeling. Liu and Hu proposed the method of influence coefficients (MIC) [23], which takes into account the elastic properties of the parts in variation analysis and springback calculation by using the pre-calculated reduced stiffness matrices. The classic MIC approach assumes that the interaction between the parts occurs only in special predefined points of fastening (e.g., welding points). This assumption is justified in many cases, but sometimes it is necessary to take into account that contact can occur at points that are not known in advance. Finding the contact points makes the analysis much more complicated because it requires solving a nonlinear contact problem. The MIC approach, combined with the direct algorithm of the contact point location, was used by Wärmefjord et al. in [24] for variation simulation in the automotive industry.
The approach based on reformulation of the contact problem in terms of quadratic programming is presented in [25]. In this paper, the variational formulation of the contact problem [26][27][28][29] is combined with implementation of a reduced stiffness matrix by using the substructuring technique (also known as Guyan reduction) [30][31][32][33]. The approach presented in [25] has been further developed and provided with practical examples in subsequent works by the same authors [22,34]. A similar approach was used by Lorin et al. [35] for simulating the assembly process in the automotive industry.
This work is devoted to the further development and extension of the numerical technique proposed in [25] for modeling and predicting the non-stationary effects that can occur during the assembly of aircraft structures. These effects include the vibration and resonance that can happen during drilling. This article also considers an example of the application of the developed technique to the simulation of drilling of the test joint and the study of accompanying vibrations with regard to variations in part shape. Figure 1 presents a sketch of a finite element model (FEM) of the assembly of two parts. An FEM consisting of a shell and bar elements is considered. The area of possible contact is called the junction area. Outside this area, there is no contact between the parts. The computational nodes in the junction area are marked with circles in Figure 1. If the possibility of contact interaction of the parts is not taken into account, the The variables of the contact problem are all degrees of freedom (DOFs) of the FEM which form a vector-valued function

Numerical Approach
where N DOF is the number of DOFs of the model and t ∈ (0, T] is the time. If the possibility of contact interaction of the parts is not taken into account, the vector-valued function U(t) obeys the equation of linear structural dynamics: where M, B and K are the matrices of mass, damping and stiffness, respectively (size N DOF × N DOF ), F(t) is the vector of the applied load (a given function of time) and the superscript dots mean differentiations in time.
The initial conditions for Equation (1) are .
where U 0 and V 0 are the given vectors of the initial displacement and velocity, respectively. The compliant parts are fixed in the assembly jig (e.g., using clamps, as in Figure 1), so some DOFs are restricted, particularly to avoid the rigid body motion of each part: where Γ f ix is the set of indices of nodes with restricted degrees of freedom. It is assumed that the finite element meshes of the parts are conformal in the junction area, so the node-to-node contact model can be used. This means that the pairs of possible contact nodes are predefined. The time-independent vector of the initial gap determines the constraints on the displacements of the parts. The non-penetration condition of Equation (4) has also been added to the system of Equations (1)-(3): where A is a time-independent matrix of size N NC × N DOF defining the constraints in the nodes. The system of Equations (1)-(4) corresponds to a contact problem with continuous time. The dimension of the problem is equal to the number of all degrees of freedom of the assembly provided by the FEM, which can reach 10 6 -10 7 variables for industrial applications. It makes this formulation of the contact problem inapplicable if serial calculations are required (e.g., variation simulation or assembly optimization) [36]. Guyan reduction [30,31] is used for the reduction of the problem dimension. This approach has been successfully applied to solving static assembly problems (see [34,37]).

Guyan Reduction
After renumbering, the vector of all DOFs can be represented as follows: where U C is the vector of those degrees of freedom in the junction area that are constrained by Equation (4). All other DOFs form the vector U R . The global stiffness matrix K can also be represented in the form of corresponding blocks: The mass matrix M, damping matrix B and load vector F can be subdivided in the same way. Then, Equation (1) can be rewritten as Usually in the assembly problems, there are no loads outside the junction area, and therefore F R = 0. As such, when excluding U R from Equation (5), the following equation for U C is obtained: It is assumed that the first and the second terms in Equation (6) can be neglected. Reduced matrices of the mass, damping and stiffness are also introduced: Then, the main equation of the original problem (Equation (1)) can be rewritten as follows: The contact problem with Equation (7) involves only constrained DOFs. Therefore, the dimension of this problem is much smaller than the original contact problem with Equation (1). Thus, this problem formulation is more suitable for the serial calculations. Note that there are various ways to reduce the mass and damping matrices, which are described in detail, for example, in [38]. Furthermore, only the reduced problem is discussed; therefore, for simplicity, the subscript C is omitted in Equation (7).

Time Discretization Algorithm
For discretization of Equation (7) in the time domain, a time grid t n = n∆t is introduced, with n = 0, . . . , N T and ∆t = T/N T and denoting U n := U(t n ). The generalized α method proposed in [39] is used. The discrete velocity V n and acceleration W n are calculated by the following formulas in accordance with the generalized α method: Here, γ and α are the parameters of the generalized α method that define the contribution of the n and (n + 1) time layers in the numerical scheme. Introducing the parameters κ and ϕ allows one to suppress the unwanted high-frequency noise during time integration. Finally, the system in Equation (10) is obtained: This can be considered a system of equations of the static contact problem with the equivalent matrix of stiffness and the equivalent load vector If the time step ∆t is constant, the matrix D does not change between steps. The load vector C n depends on the displacement, velocity and acceleration in the previous layer.
If the parameters satisfy the following conditions α γ then the generalized α method is unconditionally stable and has the second order of accuracy [39].
The system in Equation (10) with the conditions of Equations (2)-(4) for U n determines the non-linear contact problem on the time layer n. Let us formulate this problem in variational form [40] by finding the minimizer of the functional at every time layer n = 1, . . . , N T . The initial conditions for this problem are defined in Equation (2).

Basic Assumptions of the Numerical Algorithm
To summarize this section, the main assumptions of the numerical approach are listed and discussed here.
Shell and bar elements are used to build the finite element model. The fastener and drilling loads are defined as time-dependent point forces normal to the contact surface. The indicated simplifications are accepted due to the fact that the described approach is intended for modeling the assembly of large-scale aircraft structures (e.g., wing panels, fuselage sections and final aircraft assembly).
The deformation of each part in the assembly is described by the linear elasticity theory. The contact interaction between the parts without friction is taken into account. Only displacements of parts in the normal direction (to the contact surface) are modeled. A node-to-node contact model is used. These assumptions are justified, since the installed fasteners prevent significant relative tangential displacement of the parts. In addition, this work is devoted to the study of opening the gap between parts when drilling. Only normal displacements can lead to such an opening.
When modeling the non-stationary processes, the inertia and damping of the assembled parts outside the junction area are neglected. At the same time, inside the junction area, these effects are taken into account in full. This assumption makes it possible to significantly reduce the dimension of the problem by performing the Guyan reduction procedure. As a rule, the clamping of parts to the assembly stand eliminates noticeable movements (and, consequently, accelerations) of parts outside the junction area, which makes it possible to neglect the inertial effects in this area. If for some reason the inertial effects outside the junction area need to be accounted for, the developed approach cannot be applied directly and needs to be modified.
The developed approach is designed to analyze the assembly process, taking into account the deviations of the shape of the assembled parts from the nominal. Modeling of the deviations is carried out by generating the random uneven initial gaps between the parts using the stochastic methods [34]. The analysis of the influence of the part deviations is carried out by the massive solving of non-stationary contact problems and the subsequent statistical analysis of the obtained results.

Verification with Abaqus
This section is devoted to verification of the proposed approach with the commercial code Abaqus FEA (see [41] for details). A simplified assembly model was used for this. The assembly consisted of two parts: a compliant rectangular panel and a fixed absolutely rigid obstacle. The upper part was fixed to the rear edge ( Figure 2). It was modeled by the shell finite elements. A constant initial gap between parts was considered here, which equaled 6 mm everywhere in the junction area. A time-dependent harmonic load was applied to the top panel in the middle of the junction area (see Figure 2). The calculations were carried out for t from 0 to 2 s with a time step ∆t = 10 −2 s.
A comparison of the solution with Abaqus is presented for the vertical displacement of the node to which the drilling load was applied. The time interval from 0.6 s to 1 s corresponds to the time period. The displacements for this time interval are shown in Figure 3. It can be seen that the two solutions were close to each other. The slight time shift was apparently caused by inaccuracies in the reduction of the mass matrix.
Note that the upper panel came into contact with the obstacle at times t = 0.6 s and t = 1 s (the position of the obstacle is shown in Figure 3 (below) with a black line).
The computation time with Abaqus was 55.3 s, while the developed numerical procedure solved the problem in 1.2 s (excluding the time required to compute the reduced matrices). The problem was solved on a PC with an Intel Core i5 processor (3.60 GHz) and 16 GB of RAM running Windows 10. This significant speedup makes it possible to use the developed approach in variation analysis, which implies multiple computations with different input parameters.

Drilling Process Simulation
The considered model of assembly of two compliant panels is shown in Figure 4. The geometrical and mechanical parameters of the assembly were selected to imitate the part of the wing-to-fuselage joint of an aircraft. The lower panel was reinforced with two stringers (see Figure 4a). The thickness was 5 and 10 mm for the lower and upper panels, respectively. The finite element model consisted of the shell elements. The fixed edges are marked in Figure 4 with black triangles. All displacements and rotations were forbidden in all nodes of the marked edges. The panels were made of Al2024 aluminum alloy. The elastic modulus was 73 GPa, the Poisson ratio was 0.33, and the density was 2780 kg/m 3 . The junction area is depicted with green points in Figure 4. There were 231 nodes in the junction area of each part.

Drilling Process Simulation
The considered model of assembly of two compliant panels is shown in Figure 4. The geometrical and mechanical parameters of the assembly were selected to imitate the part of the wing-to-fuselage joint of an aircraft. The lower panel was reinforced with two stringers (see Figure 4a). The thickness was 5 and 10 mm for the lower and upper panels, respectively. The finite element model consisted of the shell elements. The fixed edges are marked in Figure 4 with black triangles. All displacements and rotations were forbidden in all nodes of the marked edges. The panels were made of Al2024 aluminum alloy. The elastic modulus was 73 GPa, the Poisson ratio was 0.33, and the density was 2780 kg/m 3 . The junction area is depicted with green points in Figure 4. There were 231 nodes in the junction area of each part. The calculation of the reduced matrices , and was carried out in MSC Nastran FEA code. For solving the quadratic programming problem (Equation (13)), the interior point method in MATLAB was used.
In the assembly, there were 15 holes for fastener installation (see Figure 5). Ten temporary fasteners were installed in the holes marked with yellow circles in Figure 5. As is shown in [42], the load in the temporary fasteners could be considered constant. The load in each fastener was set to 5000 N and applied to both panels as presented in Figure  5.
The panels were drilled from top to bottom. The interlayer gap was maximal when the upper panel had already been drilled and the bottom one was starting to be drilled. It was this case that was modeled in the example under consideration. The drilling load acted on the lower panel at the point marked with a red circle in Figures 4b and 5. At this point, a piecewise constant periodic load was applied with a frequency of 53.05 Hz, corresponding to the rotation speed of the drill (see the graph in Figure 4).   The calculation of the reduced matrices M C , B C and K C was carried out in MSC Nastran FEA code. For solving the quadratic programming problem (Equation (13)), the interior point method in MATLAB was used.
In the assembly, there were 15 holes for fastener installation (see Figure 5). Ten temporary fasteners were installed in the holes marked with yellow circles in Figure 5. As is shown in [42], the load in the temporary fasteners could be considered constant. The load in each fastener was set to 5000 N and applied to both panels as presented in Figure 5.
The panels were drilled from top to bottom. The interlayer gap was maximal when the upper panel had already been drilled and the bottom one was starting to be drilled. It was this case that was modeled in the example under consideration. The drilling load acted on the lower panel at the point marked with a red circle in Figures 4b and 5. At this point, a piecewise constant periodic load was applied with a frequency of 53.05 Hz, corresponding to the rotation speed of the drill (see the graph in Figure 4).    The simulation was carried out for a time interval from 0 to 0.5 s. The time step ∆t was chosen to be equal to 10 −4 s. Four random initial gaps [34] were considered to illustrate the combination of transient contact analysis and variation simulation.
The nodes of the computation mesh colored by initial gaps are plotted in Figure 6. The initial conditions for the transient simulation corresponded to the stress-strain state of the assembly with installed fastening elements. The displacements U 0 were obtained by solving the static problem with applied constant loads from the fasteners. The initial velocities V 0 were taken as equal to zero. The residual gaps between pars after installation of the fastening elements are presented in Figure 7. The simulation was carried out for a time interval from 0 to 0.5 s. The time step was chosen to be equal to 10 s. Four random initial gaps [34] were considered to illustrate the combination of transient contact analysis and variation simulation.
The nodes of the computation mesh colored by initial gaps are plotted in Figure 6. The initial conditions for the transient simulation corresponded to the stress-strain state of the assembly with installed fastening elements. The displacements were obtained by solving the static problem with applied constant loads from the fasteners. The initial velocities were taken as equal to zero. The residual gaps between pars after installation of the fastening elements are presented in Figure 7.
The aim of the simulation was to compare the oscillation amplitudes of the residual gap ( ) between the panels at the drilling point for the different initial gaps. The time dependencies of the gap are shown in Figure 8     The aim of the simulation was to compare the oscillation amplitudes of the residual gap g drill (t) between the panels at the drilling point for the different initial gaps. The time dependencies of the gap are shown in Figure 8  The values of the mean residual gap and the amplitudes of the periodic oscillations (calculated as half the difference between the maximum and minimum) are presented in Table 1 for each initial gap. It can be seen that the oscillation amplitudes depended on the initial gaps and did not correlate with the mean value of the residual gap. Consequently, non-stationary analysis can provide new information about the quality of the assembly process compared with stationary analysis.  The values of the mean residual gap and the amplitudes of the periodic oscillations (calculated as half the difference between the maximum and minimum) are presented in Table 1 for each initial gap. It can be seen that the oscillation amplitudes depended on the initial gaps and did not correlate with the mean value of the residual gap. Consequently, non-stationary analysis can provide new information about the quality of the assembly process compared with stationary analysis.

Search for Resonance Effects
This section focuses on finding possible resonance effects that can arise during the drilling process and cause dangerous vibration of the parts. This is performed by determining the natural frequencies of the assembly and testing different drilling frequencies in their neighborhood. As the stress-strain state of the assembly is controlled by the initial gap between the parts and installed fastening elements, all these factors need to be included in the analysis. Therefore, the previously described procedure for transient contact analysis was used as the main research tool.
Since it is advisable to look for resonance effects in the neighborhood of the natural frequencies, further calculations were aimed at finding them.

Determination of Natural Frequencies
The considered process was nonlinear due to the contact interaction in the assembly, and the gap between the parts could be opened and closed during drilling. Therefore, the standard technique for determination of the natural frequencies by modal analysis was not applicable here. Therefore, a different approach was used based on Fourier analysis of the system response to a sudden perturbation, which is commonly employed for examination of nonlinear systems.
Here is a brief description of the numerical algorithm:

1.
The upper panel is affected by a sudden impulse to the drilling node (indicated by the red circle in Figure 5). The impulse causes free oscillations of both panels in the assembly. The damping is set to zero, so the oscillations do not decay. These oscillations are simulated by solving Equation (13) with initial conditions corresponding to the installed fasteners, as it was described in the previous section.

2.
The frequency response function A f ree g (ω) is obtained as a Fourier transform of the computed residual gap in the drilling node. 3.
The natural frequencies ω NF of the assembly correspond to the local maxima of A f ree g (ω). The natural frequencies are calculated on a time grid with parameters ∆t = 5 × 10 −4 s and T = 1 s. These parameters allow for obtaining the natural frequencies in the interval between ω min = 10 Hz and ω max = 1000 Hz.
The resulting frequency response functions obtained by the algorithm described above are shown in Figure 9 for different initial gaps. The first five natural frequencies are marked by red rings in Figure 9 and are listed in Table 2 in ascending order. As can be seen, the calculated natural frequencies depend, although not very strongly, on the initial gap.

Building the Dependence of the Vibration Amplitude on the Drill Rotational Speed by Massive Computations
The various drilling frequencies with different initial gaps between parts and damping parameters are investigated here.
Let us define a mesh in a frequency domain that consists of 200 nodes and condenses near the natural frequencies as shown in Figure 10.

Building the Dependence of the Vibration Amplitude on the Drill Rotational Speed by Massive Computations
The various drilling frequencies with different initial gaps between parts and damping parameters are investigated here.
Let us define a mesh in a frequency domain that consists of 200 nodes and condenses near the natural frequencies as shown in Figure 10.
In this type of task, it is quite difficult to set the damping of the system correctly. Usually, this parameter is determined by experimentation on a real physical structure. The original damping matrix B was calculated with Nastran as the weighted sum of the mass and stiffness matrices [43]: where α M is the Rayleigh damping mass matrix scale factor and α K is the Rayleigh damping stiffness matrix scale factor. In order to test the influence of damping on the solution, the factors b 1 > b 2 > b 3 were introduced before the damping matrix B: b n = 10 −n , n = 1, 2, 3.
For each drilling frequency shown in Figure 10, the simulation of the drilling process was performed in the same way as described in Section 4. The four initial gaps (shown in Figure 6) and the three damping coefficients b 1,2,3 were considered. In total, 2400 transient contact simulations were performed, and the residual gap between panels g drill was computed in all cases. The resulting dependencies of the g drill oscillation amplitude on the drilling load frequency are plotted in Figure 11. As can be seen from Figure 11a,b, all high-frequency oscillations were damped with large damping coefficients. Even with small damping (Figure 11c), the large oscillation amplitudes of the residual gap were not detected in the full range of drilling load frequencies (including the neighborhood of the natural frequencies).

Building the Dependence of the Vibration Amplitude on the Drill Rotational Speed by Massive Computations
The various drilling frequencies with different initial gaps between parts and damping parameters are investigated here.
Let us define a mesh in a frequency domain that consists of 200 nodes and condenses near the natural frequencies as shown in Figure 10. In this type of task, it is quite difficult to set the damping of the system correctly. Usually, this parameter is determined by experimentation on a real physical structure. The original damping matrix was calculated with Nastran as the weighted sum of the mass and stiffness matrices [43]: where is the Rayleigh damping mass matrix scale factor and is the Rayleigh damping stiffness matrix scale factor. In order to test the influence of damping on the solution, the factors were introduced before the damping matrix : = 10 , = 1, 2, 3.
For each drilling frequency shown in Figure 10, the simulation of the drilling process was performed in the same way as described in Section 4. The four initial gaps (shown in Figure 6) and the three damping coefficients , , were considered. In total, 2400 transient contact simulations were performed, and the residual gap between panels was computed in all cases. The resulting dependencies of the oscillation amplitude on the drilling load frequency are plotted in Figure 11. As can be seen from Figure 11a,b, all high-frequency oscillations were damped with large damping coefficients. Even with small damping (Figure 11c), the large oscillation amplitudes of the residual gap were not detected in the full range of drilling load frequencies (including the neighborhood of the natural frequencies). for damping coefficient < (mm); and (c) for damping coefficient < < (mm).

Conclusions
This paper presents a methodology for determining potentially dangerous frequencies for drilling, at which a gap between the parts to be joined may open. This effect de-

Conclusions
This paper presents a methodology for determining potentially dangerous frequencies for drilling, at which a gap between the parts to be joined may open. This effect depends on the initial gap between the parts, as well as the location and properties of the fasteners.
In the examined virtual test assembly, the significant increase in the residual gap between parts was not observed in the whole considered range of drilling frequencies. This was due to the fact that in the considered example, the elastic forces dominated over the inertial ones, and noticeable resonance effects did not appear. As a rule, the resonance effects during drilling are also not revealed in industrial aircraft production. However, there is a possibility that this effect can come out in some particular cases. This should be avoided, as it can cause the injection of drilling cuttings between the contact surfaces of parts, which is considered extremely undesirable when assembling aircraft structures. Thus, it seems to be important to equip the software tools for modeling the assembly process of aircraft structures (such as ASRP [22]) with special modules to determine potentially dangerous resonance effects during drilling. Such modules can work according to the technique described in this work.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The abbreviations and symbols adopted throughout the paper are listed below: MIC Method of influence coefficients FEM Finite element model DOFs Degrees of freedom U(t) The vector of displacements and rotations in all nodes of the FEM N DOF The number of degrees of freedom t The time T The maximal time The inverse matrix of A A T The transposed matrix A M The mass matrix B The damping matrix K The stiffness matrix F(t) The vector of the applied load .

U(t)
The time derivative of U(t) ..

U(t)
The sec ond time derivative of U(t) U 0 The vector of the initial displacement V 0 The vector of the initial velocity The off-diagonal blocks of the damping matrix K CC , K RR The diagonal blocks of the stiffness matrix K CR , M RC The off-diagonal blocks of the stiffness matrix F C The vector of the load applied to constrained degrees of freedom F R The vector of the load applied to unconstrained degrees of freedom M C The reduced mass matrix B C The reduced damping matrix K C The reduced stiffness matrix ∆t The time step N T The number of time steps t n The discrete time U n The discrete displacement V n The discrete velocity W n The discrete acceleration α, γ, κ, ϕ The parameters of the generalized α method D The equivalent matrix of stiffness C n The equivalent load vector G 1 , G 2 , G 3 , G 4 The initial gaps g drill (t) The residual gap α M Rayleigh damping mass matrix scale factor α K Rayleigh damping stiffness matrix scale factor b n The damping factors