Identiﬁcation of Pre-Tightening Torque Dependent Parameters for Empirical Modeling of Bolted Joints

: The vibration characteristics of bolted structures are crucially affected by the pre-tightening torque. An approach for identifying the pre-tightening torque dependent stiffness parameters of bolted joints is proposed in this paper. Firstly, the interface of the bolted joint is characterized by the thin layer element with the isotropic material property, and the parameter value of the property is assigned relative to the distance from the center of the bolt; the inﬂuence of the bolt is ignored. Secondly, the model updating method is adopted to identify the parameters of thin layer elements using experimental data, and modal data under different values of pre-tightening torque in the range of 2 N · m~22 N · m are obtained; the torque wrench is used to determine the pre-tightening torque in the modal test. Finally, after identifying the material parameters using partial experimental data on pre-tightening torque range, the empirical equation of the interface parameters with the pre-tightening torque parameter is obtained by curve ﬁtting and the rest of the experimental data are used to verify the accuracy of the ﬁtted empirical equations. It is concluded that this method can obtain all the parameters of the equivalent thin layer elements within a certain range of pre-tightening torque, which can provide a reference for the empirical modeling of bolted structures, improve modeling efﬁciency and reﬂect the characteristic performance of real structural dynamics.


Introduction
Bolt connection has been widely used in mechanical structures; the vibration behavior of the assembly structure is closely related to the mechanical performance of the contact interface of the bolted joints [1][2][3][4]. In a well-connected joint, nonlinearity is often neglected by structural dynamics analysis due to the sufficient and appropriate pre-tightening torque. The dynamic characteristics of bolted structures always depend on the pre-tightening torque.
High-precision modeling of bolt connection is a key problem in structural design. The bolt can be well simulated by the 3D solid refinement model [5]. This modeling method can accurately analyze the deformation and stress characteristics of the inner and surrounding structures, and the effect of surface slip and contact effect on the mechanical properties of the structure can be fully considered. However, this method is too time-consuming and not conducive to engineering structures when many bolt connections exist. Under the condition that simulation precision is satisfied, some simplified models are proposed to analyze the bolt connection. Li et al. [6] proposed a bolt-spring model, which replaces the bolt model with a set of spring models. This method can improve the efficiency of computing. Vilela et al. [7] proposed a unitary model which can consider the contact interfaces, friction and preload. Fukuoka [8] utilized the three-dimensional FE model to analyze the mechanical behavior of multi-bolted joints. Yan et al. [9] used the thin-layer element method to represent the nonlinear factors in this study and verify the accuracy of thin-layer element modeling. Ahmadian et al. [10] compared the measured responses with the predictions of the model containing a parametric generic joint element. The parameters of the joint interface model are successfully identified by minimizing the difference between the measured responses and the model predictions. Yao et al. [11] proved that the material parameters of the partitioned thin-layer elements can be expressed for modeling aeroengine bolted joints. These methods can provide the deformation and stress results in good agreement with experimental results, and the effects of surface slip and contact clearance on the mechanical properties of the bolted structure can also be fully considered. However, these methods do not pay enough attention to the pre-tightening force of the bolts.
The appropriate bolt pre-tightening force will make the whole structure have better performance. The clamping conditions also have been found to affect the performance of joints. In the area of experimental study of bolt pre-tightening force, Cooper and Turvey [12] studied the fatigue test of a single bolt lap and concluded that failure loads, critical end distances and critical widths increased as the bolt clamping torque increased. Zhao et al. [13] used case studies on the equal pre-tightening force and bending moment effect to accurately predict the dynamic characteristic of a bolted assembly. In the field of numerical modeling, the effect of pre-tightening force must be taken into account for the establishment of an accurate FE model of the bolt structure. Zhao et al. [14] represented interface contact stiffness by implementing thin layer elements into the FE model and obtained the regularity of the contact stiffness in bolted interface changes with bolt preload. Ultimately, these studies establish the viewpoint that the element method is better for simulating the pre-tightening force. On this basis, the parameter identification of the finite elements of bolted structures has been widely studied. Yang and Park [15] proposed a method for inversion identification of the structural stiffness and damping of joints using frequency response function (FRF). Furthermore, in order to handle errors in test measurement, they applied this method to eliminate the noise in the original FRFs. Jiang et al. [16] identified the mechanical parameters of the contact interface, considering the uncertainty in the bolted structure by adopting the thin layer element method. Ren and Beards [17] improved the techniques for identifying joints using experimental data; these techniques are insensitive to measurement noise. Tsai and Chou [18] presented a synthesis formula to predict the FRF of the two-bolt-joint structure where the FRF was used to identify the stiffness and damping of the subject investigated. Yang et al. [19] proposed an approach for identifying the rotational stiffness and translational stiffness of the joint model by combining substructure synthesis and FRFs. Xu et al. [20] proposed a highfidelity modeling method for clamping boundary conditions. Most of these works focus on the mechanical behavior of bolted structures under fixed tightening torque. In the course of engineering a structure pre-tightening forces always change, and few studies have been carried out to analyze the vibration performance of bolted joints under varying tightening torque.
In this paper, bolt connection modeling and parameter identification in engineering are investigated. The bolt interface is modeled based on the thin layer element method. The target function and mode updating method are used to identify the thin layer element parameters. The relationships between the identified parameters and the change in pretightening torque are obtained. The resulting curves can provide guidance for the accurate modeling of the same type of bolt connections in engineering.

Thin-Layer Element
A common interface element based on the node element was proposed by Goodman, Taylor and Brekke [21]. The element formula is derived from the relative node displacement of the solid element around the interface element, as shown in Figure 1.

Thin-Layer Element
A common interface element based on the node element was proposed by Goodman, Taylor and Brekke [21]. The element formula is derived from the relative node displacement of the solid element around the interface element, as shown in Figure 1. In two-dimensional analysis, the constitutive relation or stress-relative displacement relation of the interface behavior is expressed as

Solid
where σn is the normal stress, τ is the shear stress, kn denotes the normal stiffness, ks = shear stiffness, vr and ur are the relative normal and the relative shear displacements, respectively, and [C]i is the constitutive matrix of the interface or joint element. In soil-structure interaction problems, it is usually assumed that the thickness e of the element is zero. In most problems, this formulation can provide satisfactory solutions for stick and slip modes for which the normal stress remains compressive. For some other modes such as debonding, the solutions are often unreliable. An analysis shows that the planar approximation of the element treated as a solid element can provide a satisfactory simulation of the finite-sized interface zone, and at the limit its results approach those of the zero-thickness element [22]. By defining the elements between adjacent contact bodies, simulation of the contact equivalent mechanical properties of the interface using thin layer elements can be implemented [23]. The proposed element essentially represents a solid element of small finite thickness. It is assumed that the thin layer element is generated by the solid element with eight nodes, as shown in Figure 2, and the displacement of any point of the element p(x, y, z) can be expressed as the following by the knowledge of the FE.  In two-dimensional analysis, the constitutive relation or stress-relative displacement relation of the interface behavior is expressed as where σ n is the normal stress, τ is the shear stress, k n denotes the normal stiffness, k s = shear stiffness, v r and u r are the relative normal and the relative shear displacements, respectively, and [C] i is the constitutive matrix of the interface or joint element. In soil-structure interaction problems, it is usually assumed that the thickness e of the element is zero. In most problems, this formulation can provide satisfactory solutions for stick and slip modes for which the normal stress remains compressive. For some other modes such as debonding, the solutions are often unreliable. An analysis shows that the planar approximation of the element treated as a solid element can provide a satisfactory simulation of the finite-sized interface zone, and at the limit its results approach those of the zero-thickness element [22]. By defining the elements between adjacent contact bodies, simulation of the contact equivalent mechanical properties of the interface using thin layer elements can be implemented [23]. The proposed element essentially represents a solid element of small finite thickness. It is assumed that the thin layer element is generated by the solid element with eight nodes, as shown in Figure 2, and the displacement of any point of the element p(x, y, z) can be expressed as the following by the knowledge of the FE.

Thin-Layer Element
A common interface element based on the node element was proposed by Goodman, Taylor and Brekke [21]. The element formula is derived from the relative node displacement of the solid element around the interface element, as shown in Figure 1. In two-dimensional analysis, the constitutive relation or stress-relative displacement relation of the interface behavior is expressed as where σn is the normal stress, τ is the shear stress, kn denotes the normal stiffness, ks = shear stiffness, vr and ur are the relative normal and the relative shear displacements, respectively, and [C]i is the constitutive matrix of the interface or joint element. In soil-structure interaction problems, it is usually assumed that the thickness e of the element is zero. In most problems, this formulation can provide satisfactory solutions for stick and slip modes for which the normal stress remains compressive. For some other modes such as debonding, the solutions are often unreliable. An analysis shows that the planar approximation of the element treated as a solid element can provide a satisfactory simulation of the finite-sized interface zone, and at the limit its results approach those of the zero-thickness element [22]. By defining the elements between adjacent contact bodies, simulation of the contact equivalent mechanical properties of the interface using thin layer elements can be implemented [23]. The proposed element essentially represents a solid element of small finite thickness. It is assumed that the thin layer element is generated by the solid element with eight nodes, as shown in Figure 2, and the displacement of any point of the element p(x, y, z) can be expressed as the following by the knowledge of the FE.
In the formula, x i , y i , z i , is the coordinate of the nodes, while N i is the shape function. The relationships between the element stress, the strain, and the node p are shown as follows: B and D are the geometric matrix and elastic matrix, respectively. l 1 and l 2 are the length and the width of the thin-layer elements, and e is the thickness in z-direction, which is shown in Figure 2. The thickness e is much smaller than the sizes l 1 and l 2 of the other two directions of the element. The in-plane strain components (ε x , ε y , γ xy ) are ignored under the circumstances, as well as the stress components (σ x , σ y , τ xy ) of the element. The values of partial derivatives of the form functions of thin layer elements of different sizes with respect to the local coordinates at Gaussian points are analyzed and compared in Table 1. Table 1. The value of partial derivatives of the shape functions of thin layer elements with respect to the local coordinates at Gaussian points.

Gaussian
Points As can be seen from the data in Table 1, ∂N i /∂z is greater than ∂N i /∂x and ∂N i /∂y. When the size in the z direction is much smaller than the size in the x and y directions, it can be considered that ∂Ni/∂x = ∂N i /∂y ≈ 0 and the strain component ε x = ε y = γ xy ≈ 0. In other words, only three strain components are not zero at the Gauss point. The strain component can be simplified to ε = [ε z γ yz γ zx ] T . Synthesizing the above analysis, it is assumed that the normal contact characteristics and tangential contact characteristics of the interface are independent of each other. The constitutive equation of the thin layer element which characterizes interface contact is: E n and G t are the normal elastic modulus and tangential shear modulus, respectively. If the tangential and normal contact properties are coupled, the coupling term can be added to Equation (5). If the contact properties are coupled to normals and tangents, the coupling term can be added to the constitutive relation (5). The constitutive equation is used when using isotropic constitutive material to simulate the thin layer element.
λ is the Lamé constant and G = E/2(1 + ν) is the shear modulus. As we know, the number of independent material parameters of the isotropic material is 2. When the Appl. Sci. 2021, 11, 9134 5 of 19 thickness of the cell is much smaller than the feature size of the other two directions, let ε x = ε y = ε xy ≈ 0. Finally, the constitutive equation of the material can be written as follows: In this equation, the normal elastic constant and the tangent elastic constant are not independent.
Combined with the constitutive relation from the aforementioned content, the element stiffness matrix of K can be obtained according to the principle of virtual work.
Isoparametric transformation is used to calculate the stiffness matrix K of the thin layer element. Figure 3 shows the isoparametric transformation of the thin layer element. ξ, η, and ζ are global coordinate symbols, and J is the Jacobian matrix. When the global coordinate system is consistent with the local coordinate system, the following expression is given: By using the two-node Gaussian integral method, the stiffness matrix K of the thin layer element is expressed as follows: where w ξ , w η , and w ζ are the Gaussian integral weight function. The thickness selection of thin layer elements affects the calculation results for the mechanical properties of structures. When the thickness is too large, there are too many components on the contact interface, which is inconsistent with the actual performance of the contact surface. When the thickness is too small, the Jacobian matrix tends to be a singular matrix, and the error is too large to be used to calculate the displacement strain relationship. The ratio coefficient is used as the selection condition for the thickness of thin layer element modeling, as shown below: The thickness selection of thin layer elements affects the calculation results for the mechanical properties of structures. When the thickness is too large, there are too many components on the contact interface, which is inconsistent with the actual performance of the contact surface. When the thickness is too small, the Jacobian matrix tends to be a singular matrix, and the error is too large to be used to calculate the displacement strain relationship. The ratio coefficient is used as the selection condition for the thickness of thin layer element modeling, as shown below: Desai et al. [24] proposed that when the range of R is 10-100 we can get an exact result, as studied in the context of the static contact problem where the thin layer element is based on a linear constitutive. Sharma and Desai [22] think reasonable results can be obtained when R is 5. In this paper, R in the case study is 10.
The bolted joints fix the interface; thus, the determination of the contact stiffness is the core of the dynamic analysis of the bolted structure. Sharma and Desai obtain the normal contact stiffness and tangential contact stiffness by testing the relationship between the stress and the displacement of the contact interface: where k n is the normal contact stiffness, k τ is the tangential contact stiffness, e is thickness, u r and v r are the normal displacement and tangential displacement of the thin layers, respectively, and σ and τ are normal stress and tangential stress. The interface is generated by the Hexahedron element with linear constitutive C, which is as follows: where E is the elastic Modulus, ν is the Poisson's ratio. Schmidt and Bograd [25] obtained the relationship between the equivalent shear modulus of the contact surface and the tangential stiffness of the contact surface by testing tangential force of the bolted joints' structure: where A is the real contact area and k is the tangential stiffness of the connection which correlates with the surface pressure and friction properties of the contact interface. The approximate elastic parameters of the thin layer element can be obtained by the above two methods of testing the contact stiffness. The thin layer element can provide the computation, distribution, and concentration of stresses and strain within the thin finite zones, and hence, can permit evaluation of progressive damage and failure as they occur in many engineering problems. In the case of this paper, the first four order bending modal data are used. The parameter identification method is used to determine the parameters of the thin layer element.
The key to solving parameter identification problems is optimization of the algorithm, the essence of which is an optimization problem to minimize the discrepancies between the measured and predicted parameters. The objective function and the constraint are defined as: Min p is a vector for identifying parameters. The expression of the objective function J(p) is the weighted and squared residual difference between the experimental and calculated modal parameters. The domain of the function J(p) is defined in a reasonable range of The minimum of the function is calculated in this domain. The ε is the residual difference of modal parameters, and z m , and z a (p) are the modal parameters of the test and calculation, respectively. The weighted matrix W is a diagonal matrix reflecting the relative weight of the residual difference of the modal parameters. Sensitivity analysis is used to solve the structural optimization and model updating problems with the iterative method. The j-th iteration can be described as follows: S j = W 1/2 ∂z j /∂p j represents the sensitivity matrix, which is the weighted Jacobian matrix of modal parameters. Using the numerical calculus, p j+1 , can be obtained by Equation (18). The modal parameters are taken as the objective function and the iterative optimization algorithm is used for the calculation so that the loss function is continuously reduced and the parameter p converges. Finally, the precise material parameters of the contact surface can be obtained.
Through optimization algorithm iteration, the parameters of thin layer elements can be converged and identified. The implementation procedure can be illustrated as follows: (1) Initialize j = 0, construct an initial FE model of the bolt connection using the thin-layer element; (2) Select the elastic parameters p through relative sensitivity analysis using the initial FE model; (3) Pair the experimental and numerical modal shapes using modal correlation analysis, calculating the modal assurance criterion (MAC) to achieve this; (4) Calculate the residual z m − z a j between the experimental and numerical modal data, solve the iteration format Equation (18), and obtain the variation ∆p = p j+1 − p j ; (5) If the variables ∆p are converged, go to step (6); otherwise set j = j + 1 and go to step (2), with the stop criterion p j+1 − p j ≤ 10 −6 ; (6) The exact elastic parameters are obtained.

Pre-Tightening Torque Dependent Parameters
The identification method of pre-tightening torque dependent parameters for empirical modeling of bolted joints in this paper is mainly divided into the following steps

Case Study
In this section, the effectiveness of the proposed method is verified by a bolt pretightening test and simulation analysis.

Modal Test
As shown in Figure 5, the study object is a bolted joint structure which is composed of two lap plates and four bolts. The size of the lap plate is measured at about 200 mm

Case Study
In this section, the effectiveness of the proposed method is verified by a bolt pretightening test and simulation analysis.

Modal Test
As shown in Figure 5, the study object is a bolted joint structure which is composed of two lap plates and four bolts. The size of the lap plate is measured at about 200 mm long × 80 mm wide × 8 mm deep and the connecting length is 80 mm. The lap plate is made of aluminum alloy. Bolts and nuts are made from low carbon steel of which the nominal diameter is M10. The material parameters of the above two components are shown in Table 2.   The test sample bolted structure was made by machining operation. A torque wrench was used to control the pre-tightening torque of a single bolt in the group. In this way, the modal test is carried out by hammering under different pre-tightening torque as shown in Figure 6.

Material ρ/(kg/m 3 ) E/(Pa) G/(Pa)
Aluminum alloy 2750 6.9 × 10 10 2.69 × 10 10 Low-carbon steel 7900 2.1 × 10 11 8.08 × 10 10 The test sample bolted structure was made by machining operation. A torque wrench was used to control the pre-tightening torque of a single bolt in the group. In this way, the modal test is carried out by hammering under different pre-tightening torque as shown in Figure 6.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 10 long × 80 mm wide × 8 mm deep and the connecting length is 80 mm. The lap pla made of aluminum alloy. Bolts and nuts are made from low carbon steel of which nominal diameter is M10. The material parameters of the above two components shown in Table 2.

Material ρ/(kg/m 3 ) E/(Pa) G/(Pa)
Aluminum alloy 2750 6.9 × 10 10 2.69 × 10 10 Low-carbon steel 7900 2.1 × 10 11 8.08 × 10 10 The test sample bolted structure was made by machining operation. A torque wre was used to control the pre-tightening torque of a single bolt in the group. In this way modal test is carried out by hammering under different pre-tightening torque as sh in Figure 6.   The piezoelectric accelerometers used in this modal test, had a mass of 20 g, charge sensitivity of 6.05 pC/ms 2 , and frequency range from 0.5 Hz to 5 kHz. The CL-YD-303 hammer had a reference sensitivity of 3.99 pC/N. CRAS V7.0, a vibration and dynamic signal acquisition and analysis system developed by Nanjing AnZheng Software Engineering Company, was used as the analytical instrument, and MaCras was used as the modal analysis software. To simulate free boundary conditions, a spring rope suspension was applied to counteract gravity; this method is effective in avoiding imported errors which affect the dynamic characteristics of the structure. The suspension plane is orthogonal to the test direction to avoid the effect of suspension conditions on the test results. The accelerometer is arranged at the end of the structure to avoid the vibration type node. Thirteen test points are arranged along the length direction of the lapping structure. The sampling frequency is set to 5000 Hz. The arrangement of measuring points is shown in Table 3. The pre-tightening torque T N = 2, 4, 6, 8, 10,12,14,16,18,20,22,25,28 N·m are applied to the bolt through a torque wrench. The first four order modal frequencies under different pre-tightening forces are obtained through multiple modal tests, which are all shown in Table 4. The experimental modal data are plotted and the relationship between the test frequency and the pre-tightening torque is obtained, as shown in Figure 7. As can be seen from Figure 7, with the increase of pre-tightening torque the first four order modal frequencies also increase at the same time. The modal frequencies of the first, the third and the fourth order gradually increase with the increase of the pre-tightening torque, while the second modal frequencies tend to be stable when the pre-tightening torque comes to 14 N·m.  The damping ratio of the bolted connection structure is related to many factors, su as material performance, medium condition, surface roughness, processing method, pr tightening force state and loading frequency, the influence of which is mostly reflected nonlinear relations. In addition, the measurement errors will also have an impact on t damping ratio. When measuring the frequency response functions, the measured resu of the damping ratio is larger than the actual result due to the effect of acceleromet weighting functions. The first four order damping ratios under different pre-tightenin forces were obtained, and are all shown in Table 5. Furthermore, the relationship betwe the experimental damping ratio and the pre-tightening torque was obtained as is show in Figure 8. The damping ratio of the bolted connection structure is related to many factors, such as material performance, medium condition, surface roughness, processing method, pretightening force state and loading frequency, the influence of which is mostly reflected in nonlinear relations. In addition, the measurement errors will also have an impact on the damping ratio. When measuring the frequency response functions, the measured result of the damping ratio is larger than the actual result due to the effect of accelerometer weighting functions. The first four order damping ratios under different pre-tightening forces were obtained, and are all shown in Table 5. Furthermore, the relationship between the experimental damping ratio and the pre-tightening torque was obtained as is shown in Figure 8.  As shown in Figure 8, the relationship between the damping ratio of the connection structure and the pre-tightening torque is generally non-linear, and there is no rule to be found. Meanwhile, the relationship between the first four order experimental modal shapes and the pre-tightening torque is shown in Figure 9.
As can be seen from Figure 9, the modal shapes are generally smooth. On the one hand, the connection stiffness is enough to provide good linearity to the structure; on the other hand, it is impossible to measure the detailed characteristics of the modal shapes in the connecting part due to the limited test precision. At the same time, it can be found that the mechanical properties of the bolt connection tend to be stable with the increase of the pre-tightening force.  Figure 8. Experimental damping ratio depends on the pre-tightening torque.
As shown in Figure 8, the relationship between the damping ratio of the connection structure and the pre-tightening torque is generally non-linear, and there is no rule to be found. Meanwhile, the relationship between the first four order experimental modal shapes and the pre-tightening torque is shown in Figure 9.
As can be seen from Figure 9, the modal shapes are generally smooth. On the one hand, the connection stiffness is enough to provide good linearity to the structure; on the other hand, it is impossible to measure the detailed characteristics of the modal shapes in the connecting part due to the limited test precision. At the same time, it can be found that the mechanical properties of the bolt connection tend to be stable with the increase of the pre-tightening force.

Modelling and Parameter Identification
When thin layer elements are used to model bolted structures, the role of bolts on structure is replaced by thin layer elements established on the contact surface. A schem diagram of the method is shown in Figure 10.

Modelling and Parameter Identification
When thin layer elements are used to model bolted structures, the role of bolts on the structure is replaced by thin layer elements established on the contact surface. A schematic diagram of the method is shown in Figure 10.

Modelling and Parameter Identification
When thin layer elements are used to model bolted structures, the role of bolts on the structure is replaced by thin layer elements established on the contact surface. A schematic diagram of the method is shown in Figure 10.  However, the pressure distribution on the bolt joint surface is not uniform. The pressure gradually decreases along the radial direction, as shown in Figure 11. The pressure distribution of the joint surface is shown in Equations (19)- (21).
where F b is the bolt pre-tightening force, r i is the radius of the bolt hole, r m is the maximum contact radius, r b is the radius of the bolt head, and h is the thickness of the connected piece. However, the pressure distribution on the bolt joint surface is not uniform. The pressure gradually decreases along the radial direction, as shown in Figure 11. The pressure distribution of the joint surface is shown in Equations (19)-(21).
where Fb is the bolt pre-tightening force, ri is the radius of the bolt hole, rm is the maximum contact radius, rb is the radius of the bolt head, and h is the thickness of the connected piece. We assuming that the normal dynamic stiffness Kn and tangential dynamic stiffness Kτ of the joint surface per unit area under uniform pressure are shown in Equations (22) and ( We assuming that the normal dynamic stiffness K n and tangential dynamic stiffness K τ of the joint surface per unit area under uniform pressure are shown in Equations (22) and (23) K n = α n · P (r) β n (22) where P(r) is the pressure of the joint surface, α n and β n are normal characteristic parameters of the joint, α τ and β τ are tangential characteristic parameters of the joint. The relationship between the E and G of the gradient connecting layer and K is shown in Equations (24) and (25) where A is the area of the connecting layer and h is the thickness of the connecting layer. The influence area of the bolt connection is divided into three areas, as shown in Figure 12.
Appl. Sci. 2021, 11, x FOR PEER REVIEW 16 of 21 In order to facilitate calculation and ignore the influence of the screw hole, the gradient connection layer is divided into two areas and changed into a square in this paper. The contact pressure of each area is the average pressure of each layer.

Gradient connection layer
Area 1

Area 2 Area 3
Thin layer 1 Thin layer 2 Figure 12. The influence area of the bolt connection, which is divided into three areas.
The FEM of the connection structure as shown in Figure 13 is established by ignoring the influence of the bolt weight and screw hole. The lap plate is simulated with a solid element, and the interface is simulated with an isotropic thin-layer element. The contact stiffness close to the bolt is higher than that away from the bolt, and two different isotropic constitutive relations are used to simulate the contact surface and then used to identify the parameters. The parameters to be identified near the bolt area (the red area in Figure  13) are elastic modulus E1 and shear modulus G1. For another area in the contact area, the parameters to be recognized away from the bolt are the elastic modulus E2 and shear modulus G2. The ratio coefficient of the thin-layer element is Ratio = 10, and the density of the material is 0 t/mm 3 . Figure 13. FEM of the connection structure.
It is important to note the setting of the initial value of the identification parameter. The parameters (E1, G1) of the thin layer element are of the same order of magnitude as Figure 12. The influence area of the bolt connection, which is divided into three areas.
In order to facilitate calculation and ignore the influence of the screw hole, the gradient connection layer is divided into two areas and changed into a square in this paper. The contact pressure of each area is the average pressure of each layer.
The FEM of the connection structure as shown in Figure 13 is established by ignoring the influence of the bolt weight and screw hole. The lap plate is simulated with a solid element, and the interface is simulated with an isotropic thin-layer element. The contact stiffness close to the bolt is higher than that away from the bolt, and two different isotropic constitutive relations are used to simulate the contact surface and then used to identify the parameters. The parameters to be identified near the bolt area (the red area in Figure 13) are elastic modulus E 1 and shear modulus G 1 . For another area in the contact area, the parameters to be recognized away from the bolt are the elastic modulus E 2 and shear modulus G 2 . The ratio coefficient of the thin-layer element is Ratio = 10, and the density of the material is 0 t/mm 3 . stiffness close to the bolt is higher than that away from the bolt, and two different isotropic constitutive relations are used to simulate the contact surface and then used to identify the parameters. The parameters to be identified near the bolt area (the red area in Figure  13) are elastic modulus E1 and shear modulus G1. For another area in the contact area, the parameters to be recognized away from the bolt are the elastic modulus E2 and shear modulus G2. The ratio coefficient of the thin-layer element is Ratio = 10, and the density of the material is 0 t/mm 3 . Figure 13. FEM of the connection structure.
It is important to note the setting of the initial value of the identification parameter. The parameters (E1, G1) of the thin layer element are of the same order of magnitude as the material parameters of the bolt where close to the bolt. However, the parameter selection (E2, G2) is an order of magnitude lower than the material parameters where away from the bolt. The initial value as well as the identification results are listed in Table 6. The iterative convergence curves of the identification parameters are as shown in Figure 14, when the pre-tightening torque is TN = 20 N·m. Figure 13. FEM of the connection structure.
It is important to note the setting of the initial value of the identification parameter. The parameters (E 1 , G 1 ) of the thin layer element are of the same order of magnitude as the material parameters of the bolt where close to the bolt. However, the parameter selection (E 2 , G 2 ) is an order of magnitude lower than the material parameters where away from the bolt. The initial value as well as the identification results are listed in Table 6. The iterative convergence curves of the identification parameters are as shown in Figure 14, when the pre-tightening torque is T N = 20 N·m.    Table 7 lists the modal frequency parameters before and after identification as well as the error between the above parameters and the experimental frequency. It can be seen from the identification results that the maximum updated error is no more than 2.5%. The results have comparatively higher identification accuracy.   Table 7 lists the modal frequency parameters before and after identification as well as the error between the above parameters and the experimental frequency. It can be seen from the identification results that the maximum updated error is no more than 2.5%. The results have comparatively higher identification accuracy. The elastic modulus E is greater than the shear modulus G, which reflects the characteristics of normal contact stiffness more than the tangential stiffness. The identification results (E 1 , G 1 ) near the bolt area are far greater than that of the bolt area (E 2 , G 2 ), which better reflects the actual situation. Therefore, the connection performance of the overlapping structure with multiple bolts can be well simulated by the of isotropic thin layer element.
Similarly, the material parameters of the thin layer element under different pretightening torque are obtained by the iterative solution method. Finally, the relationship between the identification parameters and the different pre-tightening torque is achieved, as shown in Figure 15. In order to verify the correctness of the method, we selected five points evenly in the interval of the fitting curve, obtained E1, G1, E2 and G2 through the parameter identification method, and put these values into the curve for verification. In Figure 16, the red dots are the verification points. It can be seen that the error between the red dots and the fitting curve is small, which proves the effectiveness of the method. In order to verify the correctness of the method, we selected five points evenly in the interval of the fitting curve, obtained E 1 , G 1 , E 2 and G 2 through the parameter identification method, and put these values into the curve for verification. In Figure 16, the red dots are the verification points. It can be seen that the error between the red dots and the fitting curve is small, which proves the effectiveness of the method. The curve obtained by this method is of high precision when a certain number of points and intervals are reasonably selected for fitting. The more points used, the more accurate the curve.
Through curve fitting, the relationship between the elastic modulus E1 and E2 and the shear modulus G1 and G2 under different pre-tightening torque was obtained. The resulting curves can guide the modeling of the same type of bolted joints. The specific operation is as follows: (i) The torque wrench can be used to determine the pre-tightening torque of the bolt. (ii) The parameter values of the thin layer element material E1, G1, E2 and G2 are drawn from the relation curves. Thus, an accurate dynamic model of the bolted joint structure could be established.

Conclusions
On the basis of the thin layer element with isotropic constitutive relation, this paper presents a method to simulate the stiffness of the joint surface of the bolt connection considering the variability of pre-tightening torque. The proposed method in this paper is verified by the pre-tightening test, and the following conclusions can be drawn: (1) Through curve fitting after the parameter identification, the relationship between the pre-tightening torque and the identification parameters can be obtained, which can guide accurate modeling of the same type of bolt connection. The curve obtained by this method is of high precision when a certain number of points and intervals are reasonably selected for fitting. The more points used, the more accurate the curve.
Through curve fitting, the relationship between the elastic modulus E 1 and E 2 and the shear modulus G 1 and G 2 under different pre-tightening torque was obtained. The resulting curves can guide the modeling of the same type of bolted joints. The specific operation is as follows: (i) The torque wrench can be used to determine the pre-tightening torque of the bolt. (ii) The parameter values of the thin layer element material E 1 , G 1 , E 2 and G 2 are drawn from the relation curves. Thus, an accurate dynamic model of the bolted joint structure could be established.

Conclusions
On the basis of the thin layer element with isotropic constitutive relation, this paper presents a method to simulate the stiffness of the joint surface of the bolt connection considering the variability of pre-tightening torque. The proposed method in this paper is verified by the pre-tightening test, and the following conclusions can be drawn: (1) Through curve fitting after the parameter identification, the relationship between the pre-tightening torque and the identification parameters can be obtained, which can guide accurate modeling of the same type of bolt connection. (2) The thin element method with linear constitutive relation makes modeling and parameter identification very efficient. It can also accurately reflect the contact performance of the joint structure of the bolt lap joint. The identification results for the contact surface material near the bolt area are far greater than for the distance away from the bolt area, which better reflects the inhomogeneity of the contact interface stiffness distribution.