An Extended Spring-mass Model of Single-lap Multi-bolt Composite Joints Considering Assembly Gap and Gap Shimming

To investigate the effect of assembly gap and shim on single-lap multi-bolt composite joint stiffness, an analytical model based on the spring-mass method was proposed, which converted the multi-bolt joints into individual single-bolt joint based on premise that there are no overlap regions of the highly stressed portions for adjacent holes. The proposed model considers the conical and spherical stress envelope and gradual elimination phenomenon of the bolt-hole clearances for multi-bolt joints. Meantime, an effective-to-equivalent gap area method was proposed to calculate the joint stiffness for situations with arbitrary assembly gap shape. Both experiment and finite element method for three-bolt joints were used to validate the proposed model with different situations of assembly gap and/or shim. The relative error of the shear stiffness between the analytical model and experiment is 0.31%, while that of the bolt stiffness is 19.8%. After that, four interested situations with different assembly gap and/or shims were discussed, and the maximum relative error of the shear stiffness between the analytical model and the finite element model is17.0%, while that of the bolt stiffness is 15.8%. Taking into account the complexity of composite material and that of assembly gap and gap shimming, the proposed analytical model is effective to predict the stiffness of the single-lap multi-bolt composite joints subjected to single-shear loading.

With the development of composite materials and its application in aerospace industry, more and more complicated composite parts are designed, which makes manufacturing accuracy guarantee more and more difficult in the assembling process. Therefore, assembly gap between the components to be assembled occurs during the assembly process because the shape distortion problems of the composite parts [1][2][3][4], which changes the load transfer and assembly stress conditions, even induces damages of the assembly components [5][6][7]. Most studies concentrate on the effect of assembly gap compensation, especially the effect of liquid shim thickness on the bearing behaviors [8][9][10][11][12][13]. However, there are few researches focusing on the assembly gap itself, like how assembly gap affect the bearing behavior and to what extent it can be. Zhai et al. [14] conducted an experimental study for single-lap single-bolt hybrid joint with considering the machined assembly gap. The results show that joint stiffness loss occurs earlier and peak load decreases as the assembly gap thickness increases. However, for multi-bolt joints, experimental method and finite element method are time-consuming and have their limitations to investigate the joint stiffness. Thus, a convenient analytical method need to be developed. Spring-mass model is one of analytical methods to predict the joint stiffness. Tate et al. [15] firstly proposed the spring-mass method to study load distribution of the joints made by isotropic materials, and it was modified by Nelson et al. [16] to account for the single-lap orthotropic composite joints. Then, McCarthy et al. introduced the bolt-hole clearance [17] of the joints under highly preload [18] into the spring-mass model, and investigated the stiffness of the multi-bolt composite joints [18]; in addition, the slope of the first linear region in load-displacement curve before bolt-hole clearance was eliminated was named as the shear stiffness, while that of the second linear region after bolt-hole clearance has been eliminated, was named as the bolt stiffness [17]. Olmedo et al. [19,20] developed an enhanced springmass model to account for the secondary bending effect of the bolt, which was based on the classic beam theory considering geometrical parameters, material properties, load path eccentricity and stacking sequence of the laminates. Yang et al. [21] extended the spring-mass model to the stiffness prediction for the single-lap single-bolt composite joints with considering assembly gap and gap shimming.
In this paper, a multi-bolt spring-mass model considering different assembly gap and/or shim, which was an extension of the Yang's model [21], was proposed to predict the shear stiffness and the bolt stiffness (named by McCarthy in reference [18]) of the single-lap multi-bolt composite joints. The conical and spherical envelope of the transverse shear stress as well as gradual elimination of the bolt-hole clearances for multi-bolt joints were considered and validated by three-bolt single-shear experiments. And effectiveto-equivalent gap area method was proposed to convert the arbitrary assembly gap shape to the basic circular shape. The finite element method (FEM) validated by the experiments was used as an assistant validation of the analytical method (AM) to discuss different assembly gap and/or shimming situations.

Analytical model for joint stiffness
In order to investigate the multi-bolt joints stiffness, a single-lap n-bolt composite joints ( 1 n  ) was selected as the research object, as shown in Fig. 1, whose geometry sizes were set referring to the ASTM D5961 [22]. In this model, basic assembly gap shape was assumed as a circular ring with assembly gap size  and assembly gap radius r  , any other assembly gap shape can be converted to the basic shape using the effective-to-equivalent gap area method (introduced later).
Multi-bolt single-lap joints with assembly gap and its spring-mass system model.
In Fig. 1 represents the maximum static friction for each bolt, which is the product of the preload Nj F and friction coefficient  . The particles can only move in the x-direction. For example, diagrams of particle 2j at the friction region and at the load take-up region (Fig. 2,) can be obtained by using the free-body diagrams and equations of motion [17,18], as given in Eq. (1) and Eq. (2).
For quasi-static in-plane loading, the accelerations   1, 2,..., 2 1 i a i n  of all the particles can be neglected so that dynamic equation can be expressed as the following equation [18], where K, x, F represent stiffness matrix, displacement vector and load vector.

Shear stiffness of the friction region
At the first linear region in load-displacement curve (friction region), stiffness matrix, displacement vector and load vector can be expressed by Eq. (4)~(6) [4], in which any particles except for the particle (2n+1) can be represented by particle (2j-1) and (2j) as long as were satisfied. The ellipsis dots in Eq. (4) represent other matrix elements, which have similar forms.
The laminate stiffness L k can be calculated by Eq. (7) referring to [18,19], where Li E is the equivalent longitudinal elasticity modulus; W and t are width and thickness of the laminate; p is the hole pitch between two adjacent holes or the distance between the hole center and the plate free end where load is applied for the holes on both sides; d is hole diameter [21].
In this study, only general situations that radius of the highly stressed portion is less than half of the   shear sh 1 sh 2 For transverse shear, elastic strain energy i U for a volume of material in laminate i may be described by Eq. (9) [18], where V , xz( ) i G and ei A are volume of the material with highly shear stressed portion, out-of-plane shear modulus and in-plane effective area.
As shown in Eq. (10), shear stiffness of the laminate shji k  is the ratio of friction force to deformation of the laminate i. In-plane effective area ei A , which satisfies conical and spherical envelope through the thickness (Fig. 4), changes as the assembly gap geometry or shimming conditions change.
(a) For single-bolt joint without assembly gap, Fig. 4 For single-bolt joint with assembly gap, Fig. 4(b), the in-plane effective area can be divided into two parts: one is the laminate region without gap in the interval      Fig. 4(d), the in-plane effective area can be divided into two parts: the first part is the laminate region without gap in the interval   In fact, actual assembly gap shape was irregular so that the effective-to-equivalent gap area method was proposed to convert the arbitrary gap shape to equivalent circular ring gap shape. The central idea is to calculate the effective gap area inside the highly stressed portion for arbitrary gap, then convert it to the equivalent gap area of circular ring shape, which is concentric with the hole. Taking the rectangle gap shape in Fig. 5 as an example (the assembly gap distributes at the right side of the vertical thick line, whose gap length is L  and gap width is b  ), the blue region means the highly stressed portion, the yellow region marked as eff S inside the real gap represents the effective gap area, while the white region inside the real gap is the ineffective gap area. And the equivalent gap area inside the dotted line is marked as   in the longitudinal direction and in the transverse direction according to laminate theory, m t is the load path eccentricity as the distance between the mid-plane of the laminates, efci L is the distance between the hole surface and the plate free end where load is applied, Fci E is the flexural modulus of the laminate [24], which can be obtained by flexural experiments [25], and ci I is the moment of inertia of the laminate.

Experiment
In order to validate the proposed analytical model, a typical single-shear experiment for three-bolt composite joints according to the ASTM D5961 was conducted. The three-bolt joints specimens were prepared with a 3-axis numerical control machine. As shown in Fig. 6 (Fig. 7). The mechanical properties of composite laminates IMS194/977-2 are given in  Fig. 7. Size of the single-lap three-bolt composite joints [22]. Table 1 The mechanical properties of the IMS194/977-2 and its equivalent engineering constants [21]. Four specimens in total were tested, experiment set-up and load-displacement curves were shown in Fig. 8. For the highly torqued single-lap multi-bolt joints, typical load-displacement curve contains three regions before damage occurs, which are friction region, transition region and load take-up region. The friction region is mainly produced by friction force between the laminates, whose stiffness is named as shear stiffness shear k . The transition region, where displacement arises without any increase in load, is determined by the bolt-hole clearances [4]. After the transition region, it comes to the load takeup region, the bolt shank begins to contact the laminates and load begin arise as the displacement increases, whose stiffness is named as bolt stiffness bolt k .

Finite element method
Considering the complexity of experimental conditions for multi-bolt composite joints with assembly gap and/or shim, finite element analysis was conducted and compared with the spring-mass model, as shown in Fig. 9. The laminate in Abaqus was modelled without the clamped region so that the length of the laminate was 4   C3D8R elements (in Abaqus) of enhanced hourglass control method were used for both laminates and fasteners. After mesh dependence studying, considering the computational cost and results accuracy, in-plane mesh size within the bolt head area was set as 0.5 mm, that of other areas was 3.0 mm, and four elements were set through the thickness for each laminate. The fasteners were all modelled as one-piece parts, which have the global element size of 0.5 mm. The laminates were assigned with homogenous mechanical properties as given in Table 1. Finite sliding and surface-to-surface discretization methods were used for interactions. The friction coefficients between the laminates were 0.42, while that between the hole walls and bolt shanks were set as 0.1, referring to the literature [28]. The friction coefficient between the laminate and bolt head has an influence on the transition region. Only when it is small enough, the slope of the load-displacement curve in the transition region can be regarded as zero. In this paper, the friction coefficients between the laminates and bolt heads were set as 0.01. In addition, soft springs with stiffness of 1 N/mm were added to limit all the displacement degrees of freedom of the fasteners, which can help improve convergence.
Two steps were set in Abaqus: (1) in the first step, the bolt loads of 5 kN (corresponding tightening torque of 7.7 N  m) were added to all fasteners, and displacement degrees of freedom (Ux, Uy, Uz) were limited for the two laminates; (2) in the second step, 'fix at current length' method (a load method in Abaqus) was used for the bolt loads, and displacement load in the x direction was applied to the reference point, which is constrained with the free side nodes of the upper laminate with kinematic coupling.

Results and discussion
Firstly, experimental results were used to validate the assumption of conical and spherical stress envelope in analytical model. The surface profiles of the three-bolt joints specimen after experiment was obtained by the digital microscope KEYENCE VHX-600E, as shown in Fig. 10. The worn areas at the outer surface (close to the bolt head) and at the contact surface (between two laminates) were marked with red circles, whose maximum diameters were measured. It can be seen that the worn areas caused by friction at the outer surface are smaller than that at the contact surface around all three holes, which indicates that the assumption of conical and spherical stress envelope in Fig. 3 is reasonable. In addition, diameter of the worn area around the hole 3 is greater than that around the hole 2, and that around the hole 1 is the smallest. Hole 3 is closest to the mobile loading end, while hole 1 is the farthest from the mobile loading end. The difference in severity of wear for three holes means that friction sliding contact times in transition region are different, which indicates that the local materials around three holes sequentially go into the transition region and bolt-hole clearances of three holes are eliminated gradually.  Secondly, single-shear experiment was conducted, as shown in Fig. 11, and the strains along the loading direction around three holes were measured by strain gauge (sensitivity coefficient 2.0±1%, resistance 120±1 Ω) and DH3816N strain measuring instrument (Donghua, China) in order to validate the gradual elimination of bolt-hole clearances for three-bolt composite joints. The strain gauges were attached to the laminate 2, which was clamped to the moving end of the testing machine. The distances between the center of the strain gauge and the center of the hole were all 11.85 mm for three holes. In Fig. 11, strain value near hole 3 is larger than that of hole 2, and strain value near hole 2 is larger than that of hole 1.
Meanwhile, it can be obviously seen that strain-time curves of hole 2 and hole 3 have gone out of the transition region, while that of hole 1 is still in the transition region, which means bolt-hole clearances of hole 1 is still not eliminated. It validates the assumption that the local materials around three holes sequentially go into the transition region and bolt-hole clearances of three holes are eliminated gradually. mm, assembly gap size was 1.0 mm) was illustrated in Fig. 14(b), for which the effective-to-equivalent gap area method must be used to convert the real gap shape to equivalent circular ring gap shape and then using Eqs. 'Shim A' situation is the shimming situation of the 'Gap A' situation, which has three individual assembly gap regions with assembly gap size of 1.0 mm, as shown in Fig. 14(c), while 'Shim B' situation compensates the assembly gap region that distributes across all three bolts (assembly gap size was 1.0 mm), as shown in Fig. 14

Conclusion
In this paper, a spring-mass stiffness model for single-lap multi-bolt composite joints considering assembly gap and gap shimming was proposed to predict the joint stiffness, including shear stiffness and bolt stiffness. The multi-bolt joints were converted into individually single-bolt joint based on the assumption that there are no overlap regions between the highly stressed portions for adjacent holes. And effective-to-equivalent gap area method was proposed to convert the arbitrary assembly gap shape to the indicates that the proposed analytical model was effective to predict the stiffness of the single-lap multibolt composite joints subjected to single-shear loading.

Availability of Data and Materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request.

Competing Interests
The authors declare that they have no competing interests.

Authors' Contributions
YY proposed the analytical model and wrote the article. BY proposed the idea of experiments. WJ assisted to conduct the simulations. DF helped in analysing the data. Figure 1 Multi-bolt single-lap joints with assembly gap and its spring-mass system model.

Figure 2
Free-body diagrams of particle 2j: (a) at friction region; (b) at load take-up region.

Figure 3
Supposed transverse shear stress envelope of the multi-bolt joints.   Rectangle assembly gap and its equivalent assembly gap shape. Damage observation and hole diameter measurement.

Figure 7
Size of the single-lap three-bolt composite joints [22].  Finite element model of three-bolt composite joints.

Figure 10
Surface pro les of the single-lap three-bolt composite joints after experiment.

Figure 11
Experimental set-up and strain-time curves of the single-lap three-bolt composite joints.

Figure 12
Model validation with experiments and FEM of the three-bolt joints without assembly gap.