**1. 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–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

**Citation:** Tian, Y.; Qian, H.; Cao, Z.; Zhang, D.; Jiang, D. Identification of Pre-Tightening Torque Dependent Parameters for Empirical Modeling of Bolted Joints. *Appl. Sci.* **2021**, *11*, 9134. https://doi.org/10.3390/app 11199134

Academic Editors: Jae Hyuk Lim and Angelo Luongo

Received: 16 August 2021 Accepted: 23 September 2021 Published: 30 September 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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.

#### **2. Basic Theory**

#### *2.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.

**Figure 1.** Interface element with relative nodal displacement.

In two-dimensional analysis, the constitutive relation or stress-relative displacement relation of the interface behavior is expressed as

$$\left\{ \begin{array}{c} \sigma\_n\\ \tau \end{array} \right\} = \left[ \begin{array}{cc} k\_n & 0\\ 0 & k\_s \end{array} \right] \left\{ \begin{array}{c} \upsilon\_r\\ \upsilon\_r \end{array} \right\} = [\mathsf{C}]\_i \left\{ \begin{array}{c} \upsilon\_r\\ \upsilon\_r \end{array} \right\} \tag{1}$$

where *σ<sup>n</sup>* 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*]*<sup>i</sup>* 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.

**Figure 2.** Thin-layer element.

$$\begin{cases} \ x = \sum\_{i=1}^{8} N\_i.x\_i \\\ y = \sum\_{i=1}^{8} N\_i.y\_i \\\ z = \sum\_{i=1}^{8} N\_i.z\_i \end{cases} \tag{2}$$

In the formula, *xi, yi, zi*, is the coordinate of the nodes, while *Ni* is the shape function. The relationships between the element stress, the strain, and the node *p* are shown as follows:

$$
\varepsilon = \mathbf{B} \cdot \mathbf{p} \tag{3}
$$

$$
\sigma = D \cdot B \cdot p \tag{4}
$$

*B* and *D* are the geometric matrix and elastic matrix, respectively.

*l*<sup>1</sup> and *l*<sup>2</sup> 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*<sup>1</sup> and *l*<sup>2</sup> 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.


As can be seen from the data in Table 1, *∂Ni*/*∂z* is greater than *∂Ni*/*∂x* and *∂Ni*/*∂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* = *∂Ni*/*∂y* ≈ 0 and the strain component *ε<sup>x</sup>* = *ε<sup>y</sup>* = *γxy* ≈ 0. In other words, only three strain components are not zero at the Gauss point. The strain component can be simplified to *ε* = [*ε<sup>z</sup> γyz γzx*] *<sup>T</sup>*. 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:

$$\left\{ \begin{array}{c} \sigma\_{\rm tr} \\ \tau\_{\rm tr} \\ \tau\_{\rm ty} \end{array} \right\} = \left[ \begin{array}{c} E\_{\rm tr} & 0 & 0 \\ 0 & G\_{\rm t} & 0 \\ 0 & 0 & G\_{\rm t} \end{array} \right] \left\{ \begin{array}{c} \varepsilon\_{\rm tr} \\ \gamma\_{\rm tr} \\ \gamma\_{\rm ty} \end{array} \right\} \tag{5}$$

*En* and *Gt* 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.

$$\begin{Bmatrix} \sigma\_{xx} \\ \sigma\_{yy} \\ \sigma\_{zz} \\ \tau\_{xy} \\ \tau\_{yz} \\ \tau\_{zx} \end{Bmatrix} = \begin{bmatrix} \lambda + 2G & \lambda & \lambda & & & \\ & \lambda + 2G & \lambda & & & \\ & & \lambda + 2G & & & \\ & & & G & & \\ & & & & G & \\ & & & & G & \\ & & & & & G \end{bmatrix} \begin{Bmatrix} \varepsilon\_{xx} \\ \varepsilon\_{yy} \\ \varepsilon\_{zz} \\ \gamma\_{xy} \\ \gamma\_{yz} \\ \gamma\_{zx} \end{Bmatrix} \tag{6}$$
 
$$G(P - \lambda P) = \epsilon\_{xy}$$

$$\lambda = \frac{G(E - 2G)}{E - 3G} = \frac{E\nu}{(1 + \nu)(1 - 2\nu)}\tag{7}$$

*λ* 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 thickness of the cell is much smaller than the feature size of the other two directions, let *ε<sup>x</sup> = ε<sup>y</sup> = εxy* ≈ 0. Finally, the constitutive equation of the material can be written as follows:

$$\left\{ \begin{array}{c} \sigma\_{z} \\ \tau\_{yz} \\ \tau\_{zx} \end{array} \right\} = \left[ \begin{array}{ccc} \lambda + 2G & 0 & 0 \\ 0 & G & 0 \\ 0 & 0 & G \end{array} \right] \left\{ \begin{array}{c} \varepsilon\_{z} \\ \gamma\_{yz} \\ \gamma\_{zx} \end{array} \right\} \tag{8}$$

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.

$$\delta \mathbf{W} = \int\_{-0}^{l\_1} \int\_{0}^{l\_2} \int\_{0}^{d} \sigma^T \delta \varepsilon \mathbf{dx} d\mathbf{y} dz = \delta \mathbf{u}\_{nodal}^T \mathbf{K} \mathbf{u}\_{nodal} \tag{9}$$

Isoparametric transformation is used to calculate the stiffness matrix *K* of the thin layer element.

$$\mathbf{K} = \int\_{V\_{\ell}} \mathbf{B}^{T} \mathbf{D} \mathbf{B} \mathrm{d}V = \int\_{-1}^{1} \int\_{-1}^{1} \int\_{-1}^{1} \mathbf{B}^{T} \mathbf{D} \mathbf{B} \mathrm{det}(\mathbf{f}) \mathrm{d}\xi \mathrm{d}\eta \,\mathrm{d}\tilde{\xi} \tag{10}$$

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:

$$\mathbf{J} = \begin{bmatrix} \partial \mathbf{x} / \partial \mathbf{\tilde{x}} & \partial y / \partial \mathbf{\tilde{x}} & \partial z / \partial \mathbf{\tilde{x}} \\ \partial \mathbf{x} / \partial \boldsymbol{\eta} & \partial y / \partial \boldsymbol{\eta} & \partial z / \partial \boldsymbol{\eta} \\ \partial \mathbf{x} / \partial \boldsymbol{\zeta} & \partial y / \partial \boldsymbol{\zeta} & \partial z / \partial \boldsymbol{\zeta} \end{bmatrix} = \begin{bmatrix} l\_1 / 2 & 0 & 0 \\ 0 & l\_2 / 2 & 0 \\ 0 & 0 & d / 2 \end{bmatrix} \tag{11}$$

By using the two-node Gaussian integral method, the stiffness matrix *K* of the thin layer element is expressed as follows:

$$\mathbf{K} = \sum\_{i=1}^{2} \sum\_{j=1}^{2} \sum\_{k=1}^{2} \mathbf{B}^{T} \left( \zeta\_{i\ast}^{x}, \eta\_{j\ast}^{z} \zeta\_{k} \right) \mathbf{C} \mathbf{B} \left( \zeta\_{i\ast}^{x}, \eta\_{j\ast}^{z} \zeta\_{k} \right) \det \left( \mathbf{J} \left( \zeta\_{i\ast}^{x}, \eta\_{j\ast} \zeta\_{k} \right) \right) w\_{\xi,i} w\_{\eta,j} w\_{\zeta,k} \tag{12}$$

where *w<sup>ξ</sup>* , *wη*, and *w<sup>ζ</sup>* are the Gaussian integral weight function.

**Figure 3.** Isoparametric transformation of the thin layer element.

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:

$$R = \frac{\max\left(l\_1, l\_2\right)}{\varepsilon} \tag{13}$$

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:

$$\begin{cases} \begin{array}{c} k\_{\tau} = \frac{c\tau}{c u\_{\overline{n}}}\\ k\_{\overline{n}} = \frac{c\sigma\_{\overline{n}}}{c v\_{\overline{n}}} \end{array} \tag{14} $$

where *kn* is the normal contact stiffness, *kτ* is the tangential contact stiffness, *e* is thickness, *ur* and *vr* 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:

$$\mathbf{C} = \frac{E}{(1+\nu)(1-2\nu)} \begin{bmatrix} 1-\nu & \nu & 0\\ \nu & 1-\nu & \nu\\ 0 & 0 & (1-2\nu)/2 \end{bmatrix} \tag{15}$$

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:

$$G = \frac{ke}{A} \tag{16}$$

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.

## *2.2. Parameter Identification*

Parameter identification methods are widely used in modeling problems [26–29]. 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:

$$\begin{cases} \text{Min } J(\mathbf{p}) = \boldsymbol{\varepsilon}^T \mathbf{W} \boldsymbol{\varepsilon} = \left\| \mathbf{W}^{1/2} (\mathbf{z}^m - \mathbf{z}^a(\mathbf{p})) \right\|\_2^2\\ \text{s.t. } \mathbf{p}\_1 \le \mathbf{p} \le \mathbf{p}\_2 \end{cases} \tag{17}$$

**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

**p**<sup>1</sup> ≤ **p** ≤ **p**2. 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:

$$\mathbf{W}^{1/2} \left( \mathbf{z}^m - \mathbf{z}\_j^a \right) = \mathbf{S}\_j \left( \mathbf{p}\_{j+1} - \mathbf{p}\_j \right) \tag{18}$$

**S***j*= **W**1/2*∂***z***j*/*∂***p***<sup>j</sup>* 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:


#### *2.3. 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


**Figure 4.** The flowchart for the parameter identification method of bolted joints under different pre-tightening torque.
