Modeling tangential friction based on contact pressure distribution for predicting dynamic responses of bolted joint structures

Abstract

Subjected to dynamic excitations, bolted joint interfaces exhibit nonlinear characteristics that significantly affect the dynamic response of assembled structures. In this paper, a novel modeling method is developed to improve the prediction accuracy of the tangential contact behavior of bolted joint interfaces. This method is based on the framework of the Iwan model and builds the relationship between the Iwan density function and the distribution of the contact pressure. It is the first time to give an explicit physical significance of the Iwan density function. The effectiveness of the proposed model is validated well by comparing with published experiments. Then the proposed model is integrated into a numerical model of a bolted joint structure and combined with the alternating frequency/time method to study its applicability in dynamics analysis. The model is fully compatible with the current state-of-the-art numerical analysis tools. The results show that the simulated interface responses are in good agreement with the experiment.

1 Introduction

Bolted joints are widely used in various mechanical structures owing to the advantages of high reliability and being easy to disassemble. Under dynamic excitations, the joint interface transmits a tangential friction force, and the joint halves undergo relative motions. The friction force exhibits nonlinear hysteresis characteristics and constitutes a closed hysteresis loop concerning interface relative displacements in case of cyclic loading. The area enclosed by the hysteresis loop represents energy dissipation induced by the tangential friction. The energy dissipation even accounts for 90% of the total structural damping [1]. Apart from the energy dissipation, the tangential friction can lead to the local stiffness softening with increasing excitation amplitudes and even bifurcation and chaos phenomena of dynamic responses [2, 3]. Therefore, the friction should be considered carefully in predicting dynamic responses of jointed structures.

The accurate prediction of the dynamic responses of bolted joint structures depends heavily on the description of the tangential friction behavior at contact surfaces. At present, many experimental studies have shown that bolted joint interfaces have the following contact behaviors [4,5,6,7,8]. For a low excitation amplitude, the force–displacement relationship is linear due to the elastic deformation of the joint surface. With increasing excitation amplitudes, the joint surface experiences a stick–slip transition and the relative sliding occurs first on the periphery of the contact area. As the excitation amplitude increases further, the sliding area gradually approaches the through hole until the entire surface slides (so-called gross slip state). Once the sliding amplitude exceeds the gap of the through hole, the bolt will be subjected to a severe shear load, and the contact stiffness will have a sudden increase. In this case, the reliability of the bolted joint is threatened.

The traditional linearized modeling method cannot accurately reproduce energy dissipation and stiffness softening characteristics of joint interfaces. Nonlinear fine-tuned finite element models can implement the prediction of the friction hysteresis behavior [9, 10], but its computational cost is high, especially in dynamic response prediction. Currently, a widely accepted approach is to establish a reliable interface constitutive model, which can accurately reproduce the friction behavior and can be easily integrated into a structural level dynamics analysis code.

Generally, there are two kinds of interface models applied to reproduce the tangential friction behavior, namely distributed element models [11,12,13,14,15] and whole interface models [7, 16,17,18,19]. The distributed element modeling approach utilizes many node-to-node contact elements to describe normal contact and tangential stick–slip friction behavior of local contact segments. All contact elements share the same contact parameters. This approach has been widely applied in the dynamic modeling of aircraft engine components [12,13,14]. Schwingshackl et al. [11] introduced this approach in dynamic analysis of bolted flange. Armand et al. [15, 20] developed a multiscale approach based on node-to-node contact elements to study the effects of roughness and fretting wear on the nonlinear dynamics of a bolted joint. However, this approach requires a large number of node-to-node elements to describe the overall interface behavior, which leads to a large number of nonlinear degrees of freedom in the dynamic models and further higher computational costs.

In comparison, the whole interface modeling approach describes the tangential friction behavior of the overall interface with only one contact model. The whole interface modeling approach is less computationally expensive in predicting the dynamic response of joint structures owing to the reduction of degrees of freedom of nonlinearity. The Iwan model [21, 22], Valanis model [23, 24], LuGre model [25] and Dahl model [26] belong to this type of model.

Due to a simple expression of the force–displacement relationship, the Iwan model has been widely used to simulate the tangential friction of bolted joint interfaces. The Iwan model was originally developed to reproduce the elastoplastic behavior of materials [21]. It is composed of numerous Jenkin’s elements [27] in series or parallel, and the sliding force of each element is determined by an assumed density function. Song et al. [19] added a linear spring to the original Iwan model to represent the residue stiffness at the gross slip stage and developed an adjusted Iwan beam element (AIBE) which can be compatible with two-dimensional linear beam finite element. Then, the nonlinear dynamics of a beam structure with bolted joints were calculated using an explicit integration formulation. Segalman et al. [16] presented a four-parameter Iwan model to reproduce an experimentally observed power-law relation between the energy dissipation per cycle and the amplitude of applied loads in case of small amplitude oscillatory loads. In this model, the force–displacement relationship depends on a truncated power-law distribution with one Dirac delta function. Based on this work, Li et al. [28] proposed a six-parameter Iwan model to consider two main phenomena: the residue stiffness and the power-law relationship of energy dissipation during the micro-slip regime. Brake [29] extended the four-parameter Iwan model to include a pinning behavior that can be observed when the applied tangential load is high. This model considers the contact between the shank of the bolt and the through hole. Biswas et al. [30] developed a two-state hysteresis model from the Iwan model to study hysteresis in bolted lap joints, which can capture minor loops formed by partial unloading. Wang et al. [31] considered the residual stiffness and the smooth transition of contact stiffness from the micro-slip to the gross slip regime and developed an improved Iwan model based on the four-parameter Iwan model. Recently, Li et al. [32, 33] developed a modified Iwan model considering the effect of normal behavior (load variation and possible intermittent separation) on the tangential friction.

In the Iwan model, the force–displacement relationship is determined by the density function of the critical sliding force of Jenkin’s elements. The Iwan-type models mentioned above assume the density function to be either the uniform distribution or the truncated power-law distribution. Although these assumptions can reproduce experimental results in some cases, why such a distribution function is used or its physical meaning was never explained.

In this work, the physical meaning of the Iwan density function is explicitly given and related to the distribution function of the contact pressure. A new interface model is then developed to simulate the tangential contact behavior at bolted joint interfaces. This model is based on the framework of the original Iwan model, and the physics-based density function derived from the distribution of the contact pressure.

This paper is organized as follows. Section 2 presents how to obtain the Iwan density function from the interface contact pressure distribution and deduces the force–displacement relationship for a cyclic loading case. Section 3 validates the proposed model by the cross-comparison between simulations and experimental results measured from a precisely controlled fretting apparatus. Section 4 applies the proposed model in a bolted joint oscillator and combines it with the alternating frequency/time domain method to predict the dynamic response of the oscillator. Section 5 highlights the accuracy and reliability of the proposed method.

2 Modeling interface friction

The proposed method starts from the contact pressure distribution of a bolted joint and converts the two-dimensional distribution to an equivalent one-dimensional distribution. Then, the one-dimensional distribution is used to calculate the Iwan density function (DF). Finally, the DF is substituted into the framework of the Iwan model to obtain the force–displacement relationship.

The Iwan model consists of a series of stick–slip models (Jenkin’s element [27]) in parallel, as shown in Fig. 1. Each Jenkin’s element is composed of a linear spring with stiffness $k_{t} /n$ and a Coulomb slider with a critical sliding force $f_{i}^{*} /n$, where $k_{t}$ is the total tangential stiffness, $n$ is the number of Jenkin’s elements and $f_{i}^{*}$ is the critical sliding force on the ith element. The critical sliding force on each element is determined by the DF $\varphi (f^{*} )$. For a monotonic loading case, the Iwan’s force–displacement relationship is

$$ T\left( \delta \right) = \int_{0}^{{k_{t} \delta }} f^{*} \varphi \left( {f^{*} } \right){\text{d}}f^{*} + k_{t} \delta \int_{{k_{t} \delta }}^{\infty } \varphi \left( {f^{*} } \right){\text{d}}f^{*} , $$

(1)

where $T(\delta )$ is the recovery force and $\delta$ the relative displacement between joint interfaces. The first integral represents all elements in the slip regime, and the second integral represents those in the stick regime [22].

Fig. 1

Scheme of the Iwan model consisting of a series of Jenkin’s elements in parallel

2.1 The distribution of contact pressure

The Iwan model is continuous, and its core, the Iwan DF, is also in the form of a continuous function. To facilitate the derivation of the DF, a continuous distribution function of the contact pressure needs to be known. However, different from sphere-on-sphere and cylinder-on-cylinder contacts, it is challenging to obtain the analytical distribution function of the contact pressure of bolted joints. In this work, a finite element (FE) model of a bolted joint with fine mesh grids is firstly used to calculate the contact pressure, and then, the simulated discrete contact pressure values are fitted by a continuous function.

Two identical plates with the dimensions $50 \times 40 \times 12 \,{\text{mm}}^{3}$ preloaded by an M8 bolt are modeled using the FE software. The nominal contact area is $1286\, {\text{mm}}^{2}$. The material is assumed to be linear elastic and the material property assigned is Young’s modulus equals 210 GPa and Poisson’s ratio equals 0.3. Each plate is meshed using 7840 8-nodes reduced integration elements (C3D8R). The interaction between contact surfaces is defined using a tangential “Penalty” friction formulation and a normal “Hard” contact. Figure 2a shows the FE model of the bolted joint. Three different preloads, 15 kN, 18 kN, and 20 kN, are applied on the bolt. Figure 2b–d shows the simulated contact pressure which decreases gradually from a maximum value on the edge of the hole to zero on the boundary.

Fig. 2

a Finite element model of an M8 bolt joint. Contact pressure distribution and relative error between finite element and approximate results: bN = 15 kN, cN = 18 kN, dN = 20 kN

The single-bolt joint has an annulus contact area surrounding the through hole. The distribution of the contact pressure is like a “cone” that the contact pressure decreases from the center to the edge, as seen in Fig. 3 where the case of 20 kN preload is set as an example.

Fig. 3

Contact pressure distribution of a single-bolt joint with a 20 kN preload

To simplify the modeling, the contact pressure distribution is approximated as a linear distribution. This approximation is widely accepted [9, 34, 35]. Figure 2 also depicts the relative errors between FE and approximate results. The relative error is defined as $\left( {p_{\text{FE}} - p_{\text{ap}} } \right)/p_{\text{ap}}$, where $p_{\text{FE}}$ represents the contact pressure obtained by the FE model, and $p_{\text{ap}}$ represents the approximate value. It can be seen that almost all relative errors are lower than $\pm 10\%$ for different preloads. Therefore, this approximation is trustworthy. Besides, the bolt preload does not influence the distribution form of the contact pressure. The approximated pressure distribution can be expressed as a function

$$ p\left( r \right) = p_{0} \left( {1 - \frac{r}{a}} \right), $$

(2)

where $p_{0}$ is the maximum pressure, $r$ the distance to contact center and $a$ the radius of the contact area. The relationship between the contact pressure and the bolt preload $N$ is

$$ N = \int_{0}^{a} p\left( r \right)2\pi r{\text{d}}r. $$

(3)

Substituting Eq. (2) into Eq. (3) yields $p_{0} = 3N/\pi a^{2}$.

The Iwan model is one-dimensional, while the contact pressure of bolted joints is distributed on a two-dimensional domain. To describe the pressure distribution (given in spatial coordinates) in the Jenkins elements coordinates, the two-dimensional pressure distribution needs to be converted to an equivalent (the bolt preload remains unchanged) one-dimensional distribution. To achieve this step, the contact pressure is integrated along the y-axis to give the normal force per contact width,

$$ p_{x} \left( x \right) = 2\int_{0}^{{\sqrt {a^{2} - x^{2} } }} p_{0} \left( {1 - \frac{{\sqrt {x^{2} + y^{2} } }}{a}} \right){\text{d}}y = p_{0} \sqrt {a^{2} - x^{2} } + \frac{{p_{0} x^{2} }}{a}{ \ln }\left( {\frac{\left| x \right|}{{a + \sqrt {a^{2} - x^{2} } }}} \right), $$

(4)

Figure 4 illustrates the dimensionality reduction of pressure distribution. In this process, the circular contact area becomes a one-dimensional contact region in the domain $- a \le x \le a$, and the maximum normal force per contact width is $p_{x0} = p_{0} a$.

Fig. 4

Dimensionality reduction of pressure distribution from two-dimension to one-dimension using an integral process

To relate the normal contact with tangential sliding friction, a concept of tangential sliding stress is proposed. The sliding stress can be regarded as a critical value of friction stress, which is defined using Coulomb friction law,

$$ \tau \left( x \right) = \mu p_{x} \left( x \right), \quad - a \le x \le a, $$

(5)

where $\mu$ is the friction coefficient.

2.2 The Iwan density function

2.2.1 Reordering the sliding stress distribution

Due to the symmetry of the sliding stress distribution, it is reordered for simplicity. As shown in Fig. 5, the contact segments having the same tangential sliding stress are stacked together. During this process, the total sliding force does not change. The distribution function of the reordered sliding stress $\tau^{*}$ is given as

$$ \tau^{*} \left( x \right) = 2\mu p_{0} \sqrt {a^{2} - x^{2} } + 2\frac{{\mu p_{0} x^{2} }}{a}{ \ln }\left( {\frac{x}{{a + \sqrt {a^{2} - x^{2} } }}} \right), \quad 0 \le x \le a. $$

(6)

Fig. 5

Reordering the sliding stress distribution and normalizing the contact width for matching the contact coordinates with the number of Jenkin’s elements

A fundamental idea of the proposed method is to build a relationship between the Iwan DF and the sliding stress distribution. For this purpose, the spatial contact region requires to correspond to the number of Jenkin’s elements, and the sliding stress distributing from the spatial domain (x) to the Jenkin’s element domain (i) is then mapped. This mapping can be implemented by normalizing the contact width, as shown in Fig. 5. The sliding stress distribution on the normalized contact width is written as

$$ \tau^{*} \left( {\frac{x}{a}} \right) = 2\mu p_{0} a\left\{ {\sqrt {1 - \left( {\frac{x}{a}} \right)^{2} } + \left( {\frac{x}{a}} \right)^{2} { \ln }\left[ {\frac{{\frac{x}{a}}}{{1 + \sqrt {1 - \left( {\frac{x}{a}} \right)^{2} } }}} \right]} \right\}, \quad 0 \le \frac{x}{a} \le 1, $$

(7)

in which $x/a$ denotes the normalized contact width, which can be replaced by $i/n$.

Then substituting $i/n$ into Eq. (7) to replace $x/a$ and multiplying by the contact width $a$, the distribution function of the critical sliding force is obtained,

$$ f^{*} \left( {\frac{i}{n}} \right) = a \cdot \tau^{*} \left( {\frac{i}{n}} \right) = 2\mu p_{0} a^{2} \left\{ {\sqrt {1 - \left( {\frac{i}{n}} \right)^{2} } + \left( {\frac{i}{n}} \right)^{2} { \ln }\left[ {\frac{{\frac{i}{n}}}{{1 + \sqrt {1 - \left( {\frac{i}{n}} \right)^{2} } }}} \right]} \right\}, \quad 0 \le \frac{i}{n} \le 1. $$

(8)

2.2.2 Definition of the density function

From a statistical viewpoint, $\varphi \left( {f^{*} } \right){\text{d}}f^{*}$ is the probability that the critical sliding force falls within the range $f^{*}$ and $f^{*} + {\text{d}}f^{*}$. Therefore, the DF $\varphi \left( {f^{*} } \right)$ of the critical sliding force can be written as

$$ \varphi \left( {f^{*} } \right) = \left| {\frac{{{\text{d}}\left( {i/n} \right)}}{{{\text{d}}f^{*} }}} \right|. $$

(9)

Due to the logarithmic term and the quadratic coupling in Eq. (8), it is difficult to get the analytical expression of the DF.

To solve this problem, we use a combination of basis functions $\left\{ {1, i/n, \left( {i/n} \right)^{2} } \right\}$ to approximate the distribution function of the critical sliding force, Eq. (8). The approximation function has the following constraints: (1) the same values at ending points with the original function, $f_{a}^{*} \left( 0 \right) = f^{*} \left( 0 \right)$ and $f_{a}^{*} \left( 1 \right) = f^{*} \left( 1 \right)$; (2) the total sliding force remains the same as the original, $\int_{0}^{1} f_{a}^{*} \left( {i/n} \right){\text{d}}\left( {i/n} \right) = \int_{0}^{1} f^{*} \left( {i/n} \right){\text{d}}\left( {i/n} \right) = \mu N$. Under the constraints of the above preconditions, the approximation function has a unique expression,

$$ f_{a}^{*} \left( {\frac{i}{n}} \right) = 2\mu p_{0} a^{2} \left[ {\left( {3 - \pi } \right)\left( {\frac{i}{n}} \right)^{2} + \left( {\pi - 4} \right)\frac{i}{n} + 1} \right]. $$

(10)

Of course, other types of or higher-order basis functions can also be employed. Figure 6 shows the dimensionless critical sliding force $f^{*} /\mu p_{0} a^{2}$ versus the normalized contact width $i/n$.

Fig. 6

Dimensionless critical sliding force versus normalized contact width

Substituting Eq. (10) and $p_{0} = 3N/\pi a^{2}$ into Eq. (9) yields the corresponding Iwan DF,

$$ \varphi \left( {f_{a}^{*} } \right) = \left| {\frac{\text{d}}{{{\text{d}}f_{a}^{*} }}\left( {\frac{i}{n}} \right)} \right| = \frac{\pi }{{6\mu N\sqrt {\frac{{\pi \left( {6 - 2\pi } \right)f_{a}^{*} }}{3\mu N} + \left( {\pi - 2} \right)^{2} } }}. $$

(11)

Figure 7 illustrates the dimensionless Iwan DF $\mu N\varphi \left( {f_{a}^{*} } \right)$ versus the dimensionless sliding force $f_{a}^{*} /\mu N$. Different from the DF of the original Iwan model, that of the proposed model is not a constant but a nonlinear function concerning the critical sliding force. In fact, if a combination of linear basis functions $\left\{ {1, i/n} \right\}$ is chosen to approximate Eq. (8), the corresponding Iwan DF is the same as that of the original Iwan model. Therefore, it can be said that the original Iwan model is a low-order form of the proposed model.

Fig. 7

Dimensionless density functions, $\mu N\varphi \left( {f^{*} } \right)$ versus dimensionless critical sliding force, $f^{*} /\mu N$

2.3 Force–displacement relationship

Once the Iwan DF is known, the relationship between the tangential friction force and the relative displacement in the case of monotonic loading can be obtained by substituting the DF into Eq. (1),

$$ T\left( \delta \right) = \frac{{\left[ {2\left( {\pi - 3} \right)k_{t} \delta - \frac{{3\left( {\pi - 2} \right)^{2} \mu N}}{\pi }} \right]\sqrt {\frac{{\pi \left( {6 - 2\pi } \right)k_{t} \delta }}{3\mu N} + \left( {\pi - 2} \right)^{2} } }}{{6\left( {\pi - 3} \right)^{2} }} + \frac{{\frac{{3\left( {\pi - 2} \right)^{3} \mu N}}{\pi } + \left( {3\pi^{2} - 21\pi + 36} \right)k_{t} \delta }}{{6\left( {\pi - 3} \right)^{2} }}. $$

(12)

The proposed model is characterized using the following parameters: the friction coefficient $\mu$, the normal load $N$, and the tangential stiffness $k_{t}$.

Subjected to cyclic excitation, the bolted joint behaves hysteresis characteristics. The corresponding force–displacement relationship forms a closed curve, which can be obtained according to Masing’s hypothesis [36],

$$ \begin{array}{*{20}l} {T_{\text{rel}} \left( \delta \right) = - T\left( {\delta_{m} } \right) + 2T\left( {\frac{{\delta_{m} + \delta }}{2}} \right),} \hfill & {\quad \dot{\delta } > 0 } \hfill \\ {T_{\text{unl}} \left( \delta \right) = T\left( {\delta_{m} } \right) - 2T\left( {\frac{{\delta_{m} - \delta }}{2}} \right),} \hfill & {\quad \dot{\delta } < 0} \hfill \\ \end{array} $$

(13)

where $\delta_{m}$ denotes the amplitude of relative displacement. Figure 8 depicts normalized hysteresis loops (the normalized tangential force $T/\mu N$ versus normalized relative displacement $k_{t} \delta /\mu N$) under different amplitudes of relative displacements. The proposed model can reproduce the transition from the stick to the slip behavior of joint interfaces.

Fig. 8

Normalized tangential force $T/\mu N$ versus normalized relative displacement $k_{t} \delta /\mu N$ under different displacement amplitudes

The area enclosed by the hysteresis loop represents energy dissipation per cycle induced by friction, which can be calculated by

$$ W = \int_{{ - \delta_{m} }}^{{\delta_{m} }} \left[ {T_{\text{rel}} \left( \delta \right) - T_{\text{unl}} \left( \delta \right)} \right]{\text{d}}\delta $$

(14)

It is well known that many experiments show that the energy dissipation per cycle during micro-slip has a power-law relationship with the excitation amplitude and corresponding power is between 2.2 and 3 [7]. Figure 9 shows the evolution of energy dissipation per cycle with increasing amplitude of normalized displacements for a cyclic loading case. In the double logarithmic coordinates, the slope of the curves is about 2.95. From this point of view, the proposed modeling method has a good performance in predicting energy dissipation. In the next section, the proposed model will be applied to solve a quasi-static problem of bolted joints and to compare with experimental results.

Fig. 9

Energy dissipation per cycle, $W$ versus amplitude of normalized displacement

3 Model validation

In this section, the proposed model is used to simulate an actual bolted joint interface and the corresponding results are compared with the experimental counterparts. Besides, the comparison among the proposed model, the original Iwan model, and the Sandia’s four-parameter Iwan model is also involved.

Eriten et al. [37, 38] developed a joint fretting apparatus for measuring friction and wear loops of mechanical joints, as shown in Fig. 10. A bolted joint with a $10 \times 17\,{\text{mm}}^{2}$ nominal rectangular contact area was studied. The joint halves are clamped by two M3 bolts and oscillated by a piezoelectric actuator (PZT). The tangential relative displacement is measured using a high-precision laser nanosensor (LNS). The tangential friction force transmitted over the joint interface is measured using a tri-axial load cell. This apparatus can implement direct force–displacement measurements that avoid the errors resulted from the postprocessing of raw data.

Fig. 10

Sketch of the testing apparatus developed by Eriten et al. [38]

Four measured hysteresis cycles under different preload (234 N, 331 N, 526 N, and 721 N) are selected to compare with simulated results. The imposed maximum tangential displacement on one of the joint halves and the friction coefficient can be found in [38]. The tangential stiffness is extracted from the measured hysteresis cycles, which is equal to the slope of the force–displacement curve in the stick regime. All contact parameters and imposed displacements for the above four cases are listed in “Appendix A”.

Figure 11 shows simulated and experimental hysteresis loops under four different bolt preloads. It can be seen that the simulated hysteresis loops by the proposed model are in good agreement with the experimental counterparts in the four cases. The difference between simulated results by the original Iwan model and the proposed model is small, which shows that the higher-order terms have less influence on the results. By comparison, the difference between simulated hysteresis loops by the four-parameter Iwan model and the experimental results is more evident. From the first two cases, we can find that in the early stages of the interface motion, the four-parameter Iwan model quickly enters the micro-slip regime, while the proposed model experiences a relatively long stick phase. While in the transition from the micro-slip regime to the gross slip regime, the proposed model takes the lead in entering the gross slip state. For the last two cases, the four-parameter Iwan model predicts a more conservative tangential friction force than the proposed model at the same imposed displacement. Comparing the shape of the simulated hysteresis curves, it can be seen that the hysteresis loop simulated the proposed model is more rigid than that by the four-parameter Iwan model.

Fig. 11

Comparison between the simulated and the experimental hysteresis loops under different normal preloads: a$N = 234\;{\text{N}}$, b$N = 331\;{\text{N}}$, c$N = 526\;{\text{N}}$, d$N = 721\;{\text{N}}$

The energy dissipation per cycle under these four preload levels are calculated, as shown in Fig. 12a. The comparison between the simulated and the experimental results shows that the proposed model has a better prediction accuracy on the energy dissipation than the four-parameter Iwan model. Figure 12b depicts the relative errors between the simulated and the experimental energy dissipation. For the first cases, the relative error of the proposed model is considerably small, lower than 1%. While for the latter three cases, the relative errors are relatively high, but still less than 10%. The relative errors of the original Iwan model are slightly higher than those of the proposed model but still acceptable. However, compared with the proposed model, the relative errors of the four-parameter Iwan model are larger for all preload levels, even more than 20% in some cases.

Fig. 12

Comparison between the simulated and the experimental energy dissipation per cycle: a energy dissipation per cycle, b relative errors of predicted energy dissipation

The good agreement of the simulated results with the experimental results validates well the effectiveness of the proposed model. Compared with the four-parameter Iwan model, the proposed model can implement a better prediction of the hysteresis loops and the energy dissipation for bolted joint interfaces. Meanwhile, the Iwan density function obtained from the contact pressure distribution is more reliable than that assumed. Therefore, it can be concluded that the proposed method is capable of simulating the tangential friction hysteresis of bolted joint interfaces with high precision.

4 Dynamic response prediction of a bolted joint structure

In this section, the proposed model is integrated into a friction oscillator system with a bolted joint to show its applicability in the dynamic analysis of jointed structures. Due to the relatively simple expression of the model, it is fully compatible with current state-of-the-art numerical analysis tools. The nonlinear dynamic response will be obtained using an alternating frequency/time domain method (AFT) [20, 39].

4.1 The bolted joint oscillator

Figure 13a shows the friction oscillator system which is made of stainless steel and consists of a frame on the left side and a lumped mass on the right side [24, 40]. The frame includes a flexible part which can ensure an expected vibration mode and a resonant frequency of the setup. There is an M10 bolt to connect the frame and the mass at the center of the setup. The nominal contact area is a $40 \times 40\, {\text{mm}}^{2}$ square surrounding the through hole. The frame and the mass are suspended horizontally on fixed support to avoid possible misalignments of the motion along the horizontal direction. A swept sine excitation is provided by an electromagnetic shaker and passes through the axle wire of the equipment. The accelerations of each part are recorded by accelerometers and used to obtain displacements by a frequency domain integral method. The transmitted interfacial friction from the left frame to the right mass can be captured from the inertia force measured on the right mass. The hysteresis loop is characterized by the relative displacement, $u_{2} - u_{B}$, and the tangential friction force on the interface, in which $u_{2}$ is the displacement of the mass and $u_{B}$ the displacement of the central part of the setup, as shown in Fig. 13a.

Fig. 13

a The bolted joint friction oscillator in [39], b the simplified three degrees of freedom model in which the bolted joint interface is simulated by the proposed models

Figure 13b illustrates a simplified model of the bolted joint oscillator, which includes three degrees of freedom. The symbols $m_{1} , m_{B} , m_{2}$ denote the lumped masses of each part, $u_{1}$ the displacement of the frame, and $k_{1}$ the stiffness of the flexible part. These parameters are listed in “Appendix B.” The friction behavior on the joint interface is simulated by the proposed contact model. The friction coefficient and the values of contact stiffness can be found in [40] and in “Appendix B,” without any parameter correction. It should be noted that the effect of residue stiffness on measured hysteresis loops is visible in this experiment. Therefore, a linear spring, $k_{0}$, is included in the contact model to represent residual stiffness. The motion equation of this friction oscillator is written as,

$$ \varvec{M}\ddot{\varvec{u}}(t) + \varvec{C}\dot{\varvec{u}}(t) + \varvec{Ku}(t) + \varvec{f}_{{{\varvec{nl}}}} \left( {\varvec{u},\dot{\varvec{u}},t} \right) = \varvec{F}_{{\varvec{e}}} (t), $$

(15)

in which $\varvec{M},\varvec{ C},\varvec{ K}$ are mass, damping, stiffness matrices, respectively, $\varvec{f}_{{{\varvec{nl}}}}$ nonlinear forces vector which depends on displacement and velocity, $\varvec{F}_{{\varvec{e}}}$ excitation forces vector. In the following simulation, since friction damping dominates structural damping in this apparatus, the material damping is ignored.

4.2 Nonlinear dynamic response

Due to the discontinuous nonlinearity of friction force, Eq. (15) cannot be solved analytically in the frequency domain. In this work, the dynamic response of the friction oscillator system is solved using the multi-harmonic balance method (MHBM) [13, 41, 42] and the AFT method. A truncated Fourier series contains only the first order of the principal vibration frequency and is used to approximate the solution, which is proved accurate enough. The details can be found in “Appendix C.”

Ref. [40] studied experimentally frequency response of this friction oscillator under different conditions. In the experiment, the bolt preload is 1.5 kN and three amplitudes of excitation, 12.5 N, 25 N, and 50 N, were considered. The simulated and experimental frequency response curves of the right lumped mass m₂ are plotted in Fig. 14. The results show typical nonlinear characteristics: the resonant frequency shifts to low frequency with increasing excitation amplitude—“stiffness softening”—and the peak of the receptance amplitude experiences a reduction. The comparison between simulated and experimental results shows that the proposed modeling method has a good prediction performance for the dynamic response of bolted joint structures.

Fig. 14

Comparison of frequency response curves between simulations and experiments for different amplitudes of excitation

Besides, two measured hysteresis loops plotted near the resonant frequency for the case of 50 N in [40] are also compared with simulated results, as shown in Fig. 15. The hysteresis loop at 336 Hz in Fig. 15a experiences a gross slip process, while another one at 340 Hz in Fig. 15b is always in the micro-slip state throughout a cycle. Overall, the simulated hysteresis loops are in good agreement with the experimental results. The visible difference between the two may be caused by the following reasons. One of the main reasons is that contact parameters may change with the reverse of loading direction, for example, the residue stiffness. The slope of simulated loops at the gross slip stage in Fig. 15a is not consistent with that of the experimental loop, which possibly caused by the instability of contact parameters due to the change of contact conditions in the process of testing. From the results, it can be seen that the slope of the measured loop at the positive gross slip state in Fig. 15a is different from that at the negative gross slip state. Another main reason is that the measured normal preload or estimated friction coefficient is not very precise, because the predicted tangential force is slightly lower than the experimental result. Other possible reasons include the effect of the redistribution of contact pressure, additional structural damping, and fretting wear.

Fig. 15

Comparison of simulated hysteresis cycles with experiments for two frequencies: a$f = 336\;{\text{Hz}}$, b$f = 340\;{\text{Hz}}$

5 Conclusions

In this paper, a contact pressure-based modeling method was presented to accurately simulate the tangential contact behavior on bolted joint interfaces. Based on the framework of the original Iwan model, this method innovatively puts forward the physical meaning of the Iwan density function and successfully relates the density function with the distribution function of contact pressure. As the authors’ knowledge, it is the first time to give an explicit physical significance of the Iwan density function. Moreover, the proposed model does not introduce new parameters, which can be characterized using conventional contact parameters: friction coefficient, tangential contact stiffness, and normal preload.

The proposed model was employed to simulate an actual bolted joint interface, and the simulated results are compared with published experimental results. The comparison shows they are in good agreement, which validates the effectiveness of the proposed method. Besides, the proposed model was also compared with Sandia’s four-parameter Iwan model. The results show that the proposed model has a better prediction of the hysteresis loops and the energy dissipation than the four-parameter Iwan model. The four-parameter Iwan model is more conservative than the proposed model when used to model bolted joints. Therefore, we can conclude that the Iwan density function obtained from the contact pressure distribution is more reliable than those prior hypothetical density functions.

The proposed model was integrated into a dynamics simulation to study its applicability in the prediction of bolted joint dynamics. The frequency response was solved using MHBM and AFT methods. The presented results show that the proposed model also accurately captures the interface nonlinear hysteresis behavior in a dynamic process. It illustrates that the proposed model can be used for the analysis of the dynamics of bolted joint structures.

Finally, it is worth noting that the proposed method can be used as a benchmark to simulate the tangential contact behavior for other contact bodies. As long as the contact pressure distribution is known, the friction hysteresis expression can be mathematically derived by this method. Furthermore, the methods also can be extended to consider the effect of fretting wear by taking the evolution of friction coefficient and the tangential contact stiffness.

Appendices

Appendix A: Contact parameters and imposed displacements in the quasi-static tests

See Table 1.

Table 1 Maximum imposed displacements and contact parameters (friction coefficient and contact stiffness) for different preload levels

Appendix B: Model parameters of the bolted joint oscillator

See Table 2.

Table 2 Model parameters of the bolted joint oscillator

Appendix C: Nonlinear dynamic response simulation

The MHBM is used to approximate dynamic response through Fourier series. Taking the displacement response, as an example, the real solution can be expanded into a series of harmonic components,

$$ \varvec{u} = \varvec{U}_{\varvec{o}} + \sum_{j = 1}^{{n_{h} }} \left[ {\varvec{U}_{\varvec{j}}^{\varvec{c}} \cos \left( {j\omega t} \right) + \varvec{U}_{\varvec{j}}^{\varvec{s}} \sin \left( {j\omega t} \right)} \right], $$

(C1)

where $\varvec{U}_{\varvec{j}}^{\varvec{c}}$ and $\varvec{U}_{\varvec{j}}^{\varvec{s}}$ are the sinusoidal and cosinoidal harmonic coefficients, $\varvec{U}_{\varvec{o}}$ the constant component of displacement individually, $j$ the order of Fourier series, $n_{h}$ the truncated order of harmonics and $\omega$ the principal vibration frequency. Similarly, the vectors of nonlinear force and excitation force can also be approximated in the same way,

$$ \varvec{F}_{\varvec{e}} = \varvec{F}_{\varvec{o}} + \sum_{j = 1}^{{n_{h} }} \left[ {\varvec{F}_{\varvec{j}}^{\varvec{c}} \cos \left( {j\omega t} \right) + \varvec{F}_{\varvec{j}}^{\varvec{s}} \sin \left( {j\omega t} \right)} \right], $$

(C2)

$$ \varvec{f}_{{\varvec{nl}}} = \varvec{f}_{\varvec{o}} + \sum_{j = 1}^{{n_{h} }} \left[ {\varvec{f}_{\varvec{j}}^{\varvec{c}} \cos \left( {j\omega t} \right) + \varvec{f}_{\varvec{j}}^{\varvec{s}} \sin \left( {j\omega t} \right)} \right]. $$

(C3)

The vibration differential equation, Eq. (15), can be transformed to an algebraic form,

$$ \varvec{Z}\left(\varvec{\omega}\right)\varvec{u} = \varvec{F}_{\varvec{e}} - \varvec{f}_{{\varvec{nl}}} , $$

(C4)

where $Z\left( \omega \right)$ is the dynamic stiffness matrix which represents the linear part of the system and can be defined as,

$$ \varvec{Z}\left(\varvec{\omega}\right) = {\mathbf{diag}}\left( {\varvec{K},\varvec{Z}_{1} , \ldots ,\varvec{Z}_{\varvec{j}} , \ldots ,\varvec{Z}_{{\varvec{n}_{\varvec{h}} }} } \right), $$

(C5)

where the $\varvec{Z}_{\varvec{j}} \varvec{ }(j = 1, 2, \ldots ,n_{h} )$ is defined as

$$ \varvec{Z}_{\varvec{j}} = \left[ {\begin{array}{*{20}c} { - \left( {\varvec{j\omega }} \right)^{2} \varvec{M} + \varvec{K}} & {\varvec{j\omega C}} \\ { - \varvec{j\omega C}} & { - \left( {\varvec{j\omega }} \right)^{2} \varvec{M} + \varvec{K}} \\ \end{array} } \right]. $$

(C6)

Equation (C6) is usually rewritten as a residual form,

$$ \varvec{R} = \varvec{Z}\left(\varvec{\omega}\right)\varvec{u} - \varvec{F}_{\varvec{e}} + \varvec{f}_{{\varvec{nl}}} , $$

(C7)

where $\varvec{R}$ is the residual. The residual equation is solved iteratively using a Newton–Raphson method for the steady response amplitude,

$$ \varvec{u}^{{\left( {\varvec{j} + 1} \right)}} = \varvec{u}^{{\left( \varvec{j} \right)}} - \left( {\left. {\frac{{\partial \varvec{R}}}{{\partial \varvec{u}}}} \right|_{{\varvec{u}^{{\left( \varvec{j} \right)}} }} } \right)^{ - 1} \varvec{R}\left( {\varvec{u}^{{\left( \varvec{j} \right)}} } \right), $$

(C8)

where $\frac{{\partial \varvec{R}}}{{\partial \varvec{u}}} = \varvec{Z}\left(\varvec{\omega}\right) + \frac{{\partial \varvec{f}_{{\varvec{nl}}} }}{{\partial \varvec{u}}}$ is the Jacobian matrix of the system.

Due to the dependent of friction force on relative displacement and velocity of contact interface and its hysteresis characteristics, it is hard to be obtained in frequency domain. Therefore, the AFT method is used to implement the mutual transformation of response signals between time and frequency domains.

After the initial response $\varvec{u}\left( \varvec{t} \right)$ in time domain is obtained, the friction force $\varvec{f}_{{\varvec{nl}}} \left( \varvec{t} \right)$ can be calculated using the proposed model. Then the friction force in time domain is transformed into frequency domain using a Fast Fourier Transform. Substituting the friction force, $\varvec{f}_{{\varvec{nl}}} \left(\varvec{\omega}\right)$, into Eq. (C4) and solving iteratively the algebraic equation to get a new steady solution, $\varvec{u}(\varvec{\omega})^{{\left( {\varvec{i} + 1} \right)}}$, the updated time domain solution, $\varvec{u}(\varvec{t})^{{\left( {\varvec{i} + 1} \right)}}$, can be obtained by an inverse Fast Fourier Transform. Subsequently, the friction force is updated according to new response. An iterative process is performed until the solution is convergent. That is, the residual $\varvec{R}$ is smaller than a given tolerance.