1 INTRODUCTION

Bolted joints have found wide spread use in many machines and structures. As basic
fastening pieces, they have direct influences on the safety and reliability of the
structural system. But in practice, delayed brittle fracture without obvious signs
occurs at the first root of the thread and the transition from the bolt head to the
shank of the bolt. This leads to serious consequences. ^{Ibrahim et al. (2005)} reviewed problems pertaining to structural dynamics
with bolted joints. It was found that the axial-load distribution was so unequal that
the serious stress concentration existed at the first thread which carried more than 30%
of the total load and the first three threads carried about 70% of the total load (^{Wang et al.,1999}, and ^{Jiang, 2006}). (^{Zhao 1994}, ^{1998}) studied the stress concentration factors (SCF)
at the roots of the thread by using the Virtual Contact Loading (VCL) method and finite
element method (FEM). ^{Fukuoka et al. (2008)} and
^{Yang et al. (2013)} used finite element method
to obtain the stress distribution in a bolt. The threads were taken as cantilever beams,
and then the axial load distribution was obtained using theory of elasticity (^{Sopwith, 1948}, and ^{Yamamoto, 1980}). Based on multi-stripe polarizer and recording
micro-densitometer, ^{Kenny et al. (1985)} and ^{Fessler et al. (1983)} used the photoelastic
stress-frozen technique to study the stress distributions both at the roots of the
thread and on the screw flanks.

The energy dissipation and the damping of bolted joints have an impact on the dynamics
of the structural system assembled by bolted joints. Therefore, it is very necessary to
study the mechanical characteristics of bolted joints. Due to wear and fretting action
of the contacting surfaces, up to 90% of the damping in structures made up of bolted
members could be supplied by the joints themselves (^{Beards, 1992}). Some researchers (^{Ibrahim et
al., 2005}, ^{Ungar, 1973}, ^{Brown, 1968}, ^{Nelson
et al., 1976}, and ^{Hanks et al., 1967})
reported that the energy dissipation should depend on the clamping pressure. Under low
pressure, the shear due to friction was small; while under high pressure, slip was
small. Therefore, there was a critical pressure so that maximum energy dissipation could
be achieved. ^{Shin et al. (1991)} and ^{Song et al. (1992)} investigated the intensity of
pressure distribution at the interfaces of bolted joints and damping characteristics.
^{Iwan et al. (1966)} and ^{Gaul et al. (2001)} used Jenkins elements in parallel to simulate a
joint. Based on the harmonic balance method, they obtained equivalent stiffness and
viscous damping for a joint. ^{Ahmadian et al.
(2007)} and ^{Jaumouillé et al. (2010)}
also used the incremental harmonic balance (IHB) method to obtain dynamic response of a
model developed for the assembled structure including the generic joint interface
element. ^{Qin et al. (2013)} investigated the
bending characteristics of the bolted disk-drum joint and calculated the steady-state
response of the jointed rotor using the harmonic balance method. ^{Luan et al. (2012)} developed a mass-spring system to obtain dynamic
response of bolted flange joints subjected to coupling vibration of transverse and
longitudinal direction. ^{Quinn (2012)} studied the
role of the mechanical joint on the overall structural response using a series-series
Iwan model. ^{Song et al. (2004)} developed modified
Iwan beam elements to simulate dynamics of beam structures with bolted joints. ^{Oldfield et al. (2005)} and ^{Ouyang et al. (2006)} used the finite element method and the
experimental method to study the dynamic behavior of a bolted joint under harmonic
loading. In addition, the hysteresis loops of the torque versus the relative angular
displacement of the joint were produced by the Masing model and the Bouc-Wen model in
their studies.

When studying dynamic characteristics of bolted joints, the stress distribution is always investigated to check whether bolts have sufficient strength. This is because once a bolt fails, the bolted joint fails and there is no need to study its dynamic characteristics. In this paper, a three-dimensional finite element model is created to simulate a bolted joint under harmonic shear displacement. In order to verify the reliability of the finite element model, the axial load distribution of the bolted joint under a preload of 38.4 kN is also studied by using an analytical method. The stress concentration factors at the roots of the thread and stress variations at ten specified points are studied. The contact condition between contacting surfaces are analyzed. The effects of the preload, the amplitude of displacement, the coefficient of friction between the two clamped plates, the loading frequency and the loading path on the dynamic response of bolted joints are investigated by means of a parametric study. In addition, the evolution of the shape of the hysteresis loop is studied. The response curve is reproduced by using the Masing model and is found to be in good agreement with that from the detailed FE model.

2 DESCRIPTION OF FINITE ELEMENT MODLE

The threads are created precisely according to the ISO standard for M12×1.75 metric
screw threads. The helix angle was so small that its influence on the axial load and
stress distribution within bolted joints could be ignored (^{Zhao, 1994}). Accounting for the symmetrical characteristics of both
model and boundary condition, only one half of the structure is modeled (see Fig. 1(a)). The bolt hole diameter of the two clamped
plates is 13 mm and there is an initial gap between the bolt shaft and the side wall of
the plates. Seven threads are created for the bolt and six threads are used in the model
for the nut. Because the first three threads carry about 70% of the total load, the mesh
in this area is the finest (see Fig. 1(b)). There
are a total of 149324 nodes and 129896 eight-node linear reduced-integration solid brick
elements in the model. Only elastic material properties are considered for all materials
involved in the bolted structure, For all materials, the modulus of elasticity,
*E*, is 210 GPa and Poisson's ratio, ν, is 0.3. The damping of the
material is not considered in this analysis.

In ABAQUS/Standard, contacting surfaces are first created, then contact pairs are defined, and last the mechanical model is created to control the behavior of contacting surfaces. In the finite element model, four contact pairs are defined: the contact between Plate 1 and the head of the bolt (Contact I), the contact between the two clamped plates (Contact II), the contact between Plate 2 and the nut (Contact III), the contact between the threads (Contact IV). The finite sliding formulation is used for Contact I and Contact II. However, considering the small-scale relative motion of the contacting bodies for Contact III and Contact IV, the infinitesimal sliding formulation is used to improve the efficiency of the analysis. The initial value of the friction coefficient, μ, is 0.15. The FE simulations conducted are listed in Table 1.

Loading case | friction coefficient μ | Preload, P0 (kN) | Displacement amplitude, A0 (mm) | Loading frequency, f (Hz) | Analytical method |

1 | 0.05 | 19.2 | 0.3 | 0 | |

2 | 0.10 | 19.2 | 0.3 | 0 | |

3 | 0.15 | 19.2 | 0.3 | 0 | |

4 | 0.20 | 19.2 | 0.3 | 0 | |

5 | 0.15 | 12.5 | 0.3 | 0 | |

6 | 0.15 | 15.0 | 0.3 | 0 | the quasi-static method |

7 | 0.15 | 25.0 | 0.3 | 0 | |

8 | 0.15 | 19.2 | 0.0025 | 0 | |

9 | 0.15 | 19.2 | 0.1 | 0 | |

10 | 0.15 | 19.2 | 0.2 | 0 | |

11 | 0.15 | 19.2 | 0.3 | 200 | |

12 | 0.15 | 19.2 | 0.3 | 1000 | the dynamic analysis method |

13 | 0.15 | 19.2 | 0.3 | 1600 | |

14 | 0.15 | 19.2 | 0.3 | 2500 |

A reference point named RP-1 is created at the center of Side I of Plate 1 (see Fig. 1 (a)). Side I is constrained with RP-1 by using coupling constraint. Coupling constraint is to make one or more degrees of freedom of two bodies have the same value. In this model, it is only used in the x direction. Both Side II and its opposite side of Plate 2 are fixed. To avoid singularity in the numerical analysis when slip occurs between the two clamped plates, the first step is to fix Side I and apply a small preload on the section that is the transition from the bolt head to the shank of the bolt. The next step is used to release Plate 1. According to the standard for bolt preload, an initial value of 19.2 kN is used in the third step. The first three steps are to create contact steadily and apply the bolt preload. To apply harmonic shear displacement on RP-1, the next five steps are created to simulate harmonic shear loading: loading in the negative direction to its maximum from zero (namely, from x = 0 to the left in Fig. 1), unload to zero, loading in the positive direction to its maximum, unload to zero; and loading in the negative direction to its maximum from zero.

After velocity reversal, the shear force at Contact II changes first quickly, and then slowly. Therefore, the initial increment size for the fifth step and the seventh step is taken to be 0.02 and that for other steps is 0.1 in the quasi-static method. The initial increment size and the maximum of increment size for all steps is one percent of the cycle in the dynamic analysis method. There is a reason why the simulation needs to be conducted under displacement control rather than load control. It is difficult to calculate the value of threshold shear force resulting in macro slip between contacting surfaces. It is therefore difficult to select a load value to simulate dynamic characteristics of the bolt. If the applied load is less than the threshold shear force, micro slip occurs between contacting surfaces, but there is no macro slip. Therefore, it cannot appropriately simulate the dynamic characteristics. Conversely, if the applied load is greater than the threshold shear force, there is macro slip. When macro slip occurs between all contact surfaces, the shear load remains constant and it will not increase with the increase of displacement. So the shear force cannot achieve the pre-set value of applied load unless the bolt gets into contact with the side wall of the holes of the two clamped plates. At this point singularity in the numerical analysis occurs and computation stops.

3 FINITE ELEMENT RESULTS AND DISCUSSION

3.1 Axial Load Distribution of Thread

To verify the finite element model (FEM), numerical results are compared with
analytical results available in the published literature. The axial load distribution
of the bolted joint was investigated by (^{Yamamoto
1970}, ^{1980}). It can be expressed by
the elastic deformation coefficients k of the bolt (index b) and nut (index n). The
axial load at any section perpendicular to the axis of the bolt or nut is given by
(^{Yang et al., 2013})

where P0 is the bolt preload, L is the thread engagement length, y is the distance from the nut bearing surface, Fy is the axial load at the section with a distance of y away from the nut bearing surface, and λ is a constant defined as

Here, Ab and An are the equivalent cross-sectional areas of bolt and nut, Eb and En are the Young modulus of bolt and nut, and β is the helix angle of threads.

Hence, the axial load of ith pitch numbered from nut bearing surface is expressed by

Here, P is the thread pitch. ηi is defined as the ratio of p(i) to P0, then

Fig.2 gives the axial-load distributions from
the analytical method (^{Yamamoto, 1980}) and the
finite element method. The finite element results show good agreement with Yamamoto's
analytical results. It is found that the first thread carried about 30% of the total
load, and the first three threads carry about 2/3 of the total load.

3.2 Stress Analysis of the Bolt

Figure 3 shows the stress distribution of the bolt under a preload of 19.2 kN. The equivalent von Mises stresses at the first three roots of the thread are found to be far beyond those at other roots of the thread of the bolt. Under harmonic shear displacement, serious damage will occur by fretting wear of the contacting surfaces of the first three threads. The maximum stress value is 1091 MPa and it occurs at the first root of the thread. It is much larger than the yield stress of carbon steel.

The stress concentration factor (SCF) is defined as the ratio of the principal stress at the bottom of the root of a thread to the nominal tensile stress on the section of the bolt. The distribution of the SCF along the bolt threads is shown in Figure 4. The first root of the thread has the greatest stress concentration and hence firstly undergoes the plastic deformation. The SCF in the next roots reduce first quickly, then slowly. The SCF reduces quickly at the seventh root of the thread, this is because the last thread is a free end, not bearing bolt pretension. The stress at the first root of the thread exceeds the yield stress of bolt material. Even if the stress does not exceed the yield limit, under harmonic loading, it is still prone to fatigue fracture.

In order to study the variation of the equivalent von Mises stress with time, one half of the first thread of the bolt shown in Fig.5 in to be studied in detail. Two paths, denoted by a and b, are defined in the figure. Path a is the bottom of the root of the first thread, and Path b is the edge of the contact zone for the root of the first thread. At the same time, five planes are created perpendicular to the circumferential direction of the bolt, and the angle between every two adjoining planes is 45 degrees. Consequently, there are ten intersection points (one on path a and one on path b on each plane), denoted by a1, a2, ..., a5 and b1, b2, ..., b5. The direction of the harmonic displacement is parallel to the x-axis.

Figure 6 a) shows the stress distributions along Path a at six time moments during one and 1/4 loading cycles. After the application of preload (step 3), the transverse displacement is zero and the von Mises stresses along Path a are in uniform distribution. While the transverse displacement returns from the maximum (x = 0.3 mm) or minimum (x = ‒0.3 mm) to zero (step 5 or step 7), the stresses are not in uniform distribution. This illustrates that hysteresis phenomena maybe occur under harmonic shear displacement. When the transverse displacement reaches its maximum (step 4), the von Mises stress at Node a1 is increased least, while that at Node a5 is reduced least. Conversely, when the transverse displacement reaches its minimum (step 6), the von Mises stress at Node a1 is reduced slightest and that at Node a5 is increased slightest. The variations of the von Mises stress at the five nodes during the application of the harmonic shear displacement are shown in Figure 6 b). The stress amplitudes of Nodes a1 and a5 are almost identical with a value of 1600 MPa and larger than those of the other three nodes. Therefore, fatigue fracture would happen at Node a1 or Node a5. When the transverse displacement reaches ‒0.30mm from ‒0.18 mm, the von Mises stresses at the five nodes remain constant, which indicated the shear load is equal to the friction force and that the bolt is in balance. Figure 7 shows the stress distributions along Path b and the variations of the stress at the five nodes during the first one and 1/4 loading cycles. Compared with figure 6, similar conclusions can be obtained. The stress amplitudes of Nodes b1 and b5 are larger than those of other three nodes, but smaller than those of Nodes a1 and a5. Accordingly, under harmonic displacement, fatigue fracture is most likely to happen at the bottom of the first thread.

3.3 Contact State Analysis

^{Mindlin (1949)} reported that there were stick
and slip regions when micro slip occurred between contacting surfaces. The elastic
deformations of two balls compressed towards each other in the presence of a
tangential force and a torsional couple were first described in his work. In order to
analyze the stick-slip conditions for Contact I and Contact II, function f, the ratio
of the friction stress to the shear stress, is defined by

where μ is the friction coefficient between contacting surfaces, σ is the normal stress and τ is the shear stress on the contact surface.

When f is larger than 1, this node experiences sticking. When f is equal to 1, this node undergoes slipping. Theoretically, it is impossible for f to become smaller than 1. Figure 8 shows the variations of the maximum and minimum of f (max(f) and min(f)) for Contact I and Contact II during the loading process. When the transverse displacement returns from its maximum to zero (step 7), max(f) is greater than 1 and min(f) is equal to 1, which indicated that Contact I is in the partial slip condition.

When the transverse displacement is reducing from ‒0.18 mm to ‒0.30 mm, max(f) and min(f) are equal to 1. It is clear that Contact I is under the gross slip condition, the bolt is in balance and the shear load remains constant (Fig.8 a)). Similarly, when the transverse displacement reaches its minimum (step 6), max(f) and min(f) are equal to 1, which indicated that Contact II is under the gross slip condition. After velocity reversal, Contact II gets into the partial slip state. When the transverse displacement is reducing from 0.30 mm to 0.279 mm, max(f) and min(f) are equal to 1, which indicated that Contact II is under the gross slip condition (Fig. 8 b).

3.4 The Hysteresis Loop

A bolted joint displays a hysteresis loop for the transverse load versus the relative displacement of the joint when subjected to harmonic shear displacement. Figure 9 shows the hysteresis loops for different friction coefficients between the two clamped plates. The energy dissipation and the maximum of the shear load increase with the increase of friction coefficient. When the relative transverse displacement between the clamped plates is zero, the shear load is not zero. To explain this hysteresis phenomenon, the effects of friction coefficient for Contact I and Contact II (μb and μp) on the hysteretic curves are analysed (Figure 10). Curve I is the hysteresis loop of the joint under the condition of neglecting the friction coefficient between the two clamped plates. If the absolute value of |x - xrev| is greater than or equal to 0.48 mm, Contact I is in the gross slip condition, and the shear load ( ) remains constant. If not, Contact I is in the partial slip condition. Curve II is the hysteresis loop of the joint under the condition of neglecting the friction coefficient between Plate 1 and the head of the bolt. Similarly, if the absolute value of |x - xrev| is greater than or equal to 0.021mm, Contact II is in the gross slip condition, and the shear load ( ) is a constant. If not, Contact II is in the partial slip condition. At the moment of velocity reversal, the contact states for Contact I and Contact II transform from the gross slip condition to the partial slip condition. Contact II goes into the gross slip condition before Contact I. After the displacement reversal, the rate of change of the shear load denoted by is small, but that of the shear load named is large. The relationship between and is as follows:

Figure 11 shows the hysteresis loops for
different bolt preloads. When a preload of 12.5 kN is used, the friction forces
between contacting surfaces are small. If |*x* - *x*
*
_{rev}
*| ≥ 0.33 mm, Contact I and Contact II are under the gross slip condition,
and the shear load remains constant. With the increasing preload, the friction forces
increase and the shear load leading to slip between the plates also increases. At the
same time, the energy dissipation increases. When the preload reaches 25 kN, under
harmonic shear displacement, Contact I stays in the partial slip condition. The shape
of the hysteresis loop transforms from a parallel hexagon to a parallelogram.

Figure 12 shows the hysteresis loops for different amplitudes of the harmonic shear displacement applied to the bolt. When the displacement amplitude is 2.5 μm, Contact I and Contact II stays in the partial slip condition and the hysteresis loop appears linear. With the displacement amplitude increases, the gross slipping condition begins to occur between the two clamped plates (Contact II) and the hysteresis loop appears like a parallelogram. The energy dissipation increases at the same time. When the displacement amplitude is 0.3mm, if |x - xrev| ≥ 0.48 mm, then Contact I is under the gross slip condition. The hysteresis loop appears as a parallel hexagon and the energy dissipation is further increased.

The above results are calculated with the quasi-static method. At low frequency, the inertial force is so small that it can be ignored, and the quasi-static method can be used to study the dynamic behavior of a bolted joint under harmonic loading; while at high frequency, the inertial force is too large to ignore its effect on the hysteresis loop, and the dynamic analysis method should be used. In order to study the effect of frequency on the hysteresis loop, sinusoidal displacement is applied on RP-1. Figure 13 shows the hysteresis loops for different frequencies of the harmonic shear displacement applied to the joint. When the frequency is 200 Hz, the maximum value of the inertial force is 0.037 kN. The shear load is mainly affected by the friction force and the hysteresis loop coincides with that calculated with the quasi-static method. The effect of the inertial force on the hysteresis loop increases with the increase of frequency.

Figure 14 shows the hysteresis loops for different loading paths. At low frequency, the hysteresis loop for triangular wave loading shows good agreement with that for sinusoid wave loading; while at high frequency, the hysteresis loops for the previously mentioned loading paths are different. Consequently, at low frequency, the effect of loading path on the hysteresis loop is small; while at high frequency, the effect is remarkable.

Please note that it takes a long time to complete a single run of dynamic analysis using ABAQUS/Standard that produces the hysteresis loops in Figs. 13 and 14. This suggests that a simplified model of the bolted joint should be very useful.

4 MASING MODEL

A Jenkins element is a spring and a Coulomb element connected in series. With a single
Jenkins element the joint only has two physical states, sticking and total slipping.
Therefore, a single Jenkins element cannot appropriately simulate the state of motion of
bolted joints. In order to simulate micro-slip, Gaul and co-workers (2001) made good use
of Masing model to simulate a joint. Masing model was further developed to describe the
elasto-plastic behavior of metals (^{Masing, 1923}).
When some of Jenkins elements are sticking and others are slipping, the model is in the
partial slip condition. While all the Jenkins elements are slipping, the model is in the
gross slipping condition (Fig. 15).

The total joint force is given by (adapted from Ref. (^{Gual, 2001}))

The formulation represents a non-smooth function, so the dynamic model is nonlinear, where

and ki (i = 0, 1, 2, ..., n) is the spring constant, Ci (i = 1, 2, 3, ..., n) is the threshold force of Coulomb element 'i', x is the displacement of a Jenkins element, xrev is the displacement immediately prior to velocity reversal. The dot over the dependent variable represents the derivative with respect to time. Function sgn ( ) is defined as

The parameters k0, ki, and Ci in the Masing model are used to fit the hysteresis curve.
If selecting an nth order model, 2n + 1 equations are needed to solve these parameters.
Using Masing model, (^{Oldfield et al. 2005}) and
(^{Ouyang et al. 2006}) obtained the hysteresis
loops of the torque versus the relative angular displacement of a bolt joint. In
addition, the method for determining parameters k0, ki, and Ci was discussed in their
investigation.

The solution procedure is as follows:

1. A hysteresis loop is firstly generated by the detailed finite element model. The loop, from one point of velocity reversal to another, is broken down into n + 1 sections of apparently straight lines. At the point of the velocity reversal at the bottom of the loop, it is considered that no Jenkins element has reached its threshold force. Namely, all of the Jenkins elements are in sticking state. At each subsequent section, the contact state for a single Jenkins element transforms from sticking state to slipping state. At the last section, all of the Jenkins elements are in slipping state. Consequently, the following equations are obtained from Eq. (7) as

2. By calculating the gradient of each segment of the curve, the cumulative stiffness of all the springs active in the system is found.

where fn is the functional expression of the nth (n = 1, 2...) segment of the loop and grad fn is the gradient of that segment.

3. Substituting the known spring constant into equations (10), the resistance and threshold forces of Coulomb elements are found.

Figure 16 shows the hysteresis loops generated in the Masing model. One Jenkins element only provides a very coarse representation. However, selecting a fourth order model, the hysteresis loop almost overlaps with that predicted by the finite element method. Therefore, the response curve can be reproduced by using the fourth order Masing model.

5 CONCLUSIONS

The dynamic behavior of a bolted joint subjected to harmonic shear displacement is studied in this paper. The following conclusions can be made from the numerical results.

1. Under the bolt preload, the equivalent von Mises stress and the stress concentration factor at the first root of the thread are found to be largest, and there is a risk of fatigue fracture.

2. For the contact between the top plate and the head of the bolt, the gross slip condition will occur during a loading cycle when the preload is small and the harmonic shear displacement amplitude is large. So do the contact between the two clamped plates. After velocity reversal, the contact conditions for the two contact pairs transform from the gross slip state to the partial slip state. Compared with the contact between the top plate and the head of the bolt, gross slipping between the two clamped plates will take place earlier. The shear load remains constant after the two contact pairs undergo gross slipping.

3. The threshold force leading to slip between the two clamped plates and the energy dissipation increase with the increase of friction coefficient between the clamped plates. Similarly, they become larger when preload increases. When preload reaches a certain value, the contact between the top plate and the head of the bolt is always in the partial slip state, and the shape of the hysteresis loop transforms from a parallel hexagon to a parallelogram. When shear displacement amplitude is so small that the two contact pairs described earlier are always in the partial slip state, the hysteresis loop appears as a straight line. With increasing displacement amplitude, gross slipping firstly begins to occur between the two clamped plates, and then between the top plate and the head of the bolt. The hysteresis loop appears as a parallelogram, and then a parallel hexagon. Under low frequency, the effects of frequency and loading path are small, while under high frequency, they are remarkable.

4. The hysteresis loops of the bolted joint can be well reproduced by using a fourth order Masing model.

Nomenclature

, = transverse load calculated with quasi-static method and dynamic analysis method

*f* = frequency of the harmonic shear displacement

*k*
_{0} = permanent spring constant

*k*
*
_{i }
*= constant of the spring attached to Coulomb element '

*i*'

*x*, =
displacement and velocity of Jenkins element

*x*
*
_{rev }
*= displacement immediately prior to velocity reversal

*C*
*
_{i}
* = threshold force of Coulomb element '

*i*'

*y* = distance from the nut bearing surface

*f*
*
_{y}
* = axial load at the section with a distance of y away from the nut bearing
surface

*P*
_{0} = initial clamping force or preload

*L* = thread engagement length

*A*
*
_{b}
*,

*A*

*= equivalent cross-sectional areas of bolt and nut*

_{n}
*E*
*
_{b}
*,

*E*

*= Young modulus of bolt and nut*

_{n}β = helix angle of threads

*p*(*i*) = axial load of nth pitch numbered from nut
bearing surface

*P* = thread pitch

σ = normal pressure on the contact surface

τ = shear stress on the contact surface

μ*
_{p}
* = friction coefficient between the two clamped plates

μ*
_{b}
* = friction coefficient between the moving plate and the head of the bolt