Numerical and Theoretical Studies of Bolted Joints Under Harmonic Shear Displacement

A three-dimensional finite element model used to simulate the bolted joint is created using ABAQUS package. The stress concentration factors at the roots of the thread are first studied with a preload of 38.4 kN. Under harmonic transverse shear displacement, not only the stress variations at ten specified points of the first thread but also contact conditions between the contacting surfaces are studied. By changing the preload value, relative displacement, the coefficient of friction between the two clamped plates successively and the loading frequency, their effects on the hysteresis loops of the transverse load versus the relative displacement of the joint are then analyzed. Finally the hysteresis loops are produced by the Masing model. It is found that due to the greatest stress concentration factor, fatigue failure will occur at the root of the first thread. With the increase of preload value, the amplitude of displacement and the coefficient of friction, frictional energy dissipation of the joint increases. Very good agreement is achieved between the hysteresis loops produced from the fourth-order Masing model and the detailed ABAQUS model and therefore the fourth-order Masing model can replace the time-consuming finite element model.


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, andJiang, 2006).Zhao (1994Zhao ( , 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, andYamamoto, 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.
Latin American Journal of Solids and Structures 12 (2015) 115-132 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 threedimensional 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.

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

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 (1970Yamamoto ( , 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 P 0 is the bolt preload, L is the thread engagement length, y is the distance from the nut bearing surface, F y is the axial load at the section with a distance of y away from the nut bearing surface, and is a constant defined as Latin American Journal of Solids and Structures 12 (2015) 115-132 Here, A b and A n are the equivalent cross-sectional areas of bolt and nut, E b and E n 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   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.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.), 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 a 1 is increased least, while that at Node a 5 is reduced least.Conversely, when the transverse displacement reaches its minimum (step 6), the von Mises stress at Node a 1 is reduced slightest and that at Node a 5 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 a 1 and a 5 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 a 1 or Node a 5 .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 b 1 and b 5 are larger than those of other three nodes, but smaller than those of Nodes a 1 and a 5 .Accordingly, under harmonic displacement, fatigue fracture is most likely to happen at the bottom of the first thread.

Stress Analysis of the Bolt
Latin American Journal of Solids and Structures 12 (2015) 115-132  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

Contact State Analysis
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.

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   sing 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 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 rev -0.48 xx  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 ABA-QUS/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.

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 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 k i (i=0,1,2, …,n) is the spring constant, C i (i=1,2,3,…,n) is the threshold force of Coulomb element 'i', x is the displacement of a Jenkins element, x rev is the displacement immediately prior to velocity reversal.The dot over the dependent variable represents the derivative with respect to time.Function

 
sgn x is defined as The parameters k 0 , k i , and C i 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 k 0 , k i , and C i 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 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 f n is the functional expression of the nth (n=1, 2…) segment of the loop and grad f n 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.

Figure 2 :
Figure 2: Comparison of axial-load distributions from the analytical method and the FEM.

Figure 3
Figure3shows 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

Figure 3 :
Figure 3: Stress distribution on the threads.

Figure 4 :
Figure 4: The SCF distribution along the roots of the thread.

Figure 5 :
Figure 5: Definition of the two paths on the first thread.

Figure 6
Figure6a) 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.3mm) 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 a 1 is increased least, while that at Node a 5 is reduced least.Conversely, when the transverse displacement reaches its minimum (step 6), the von Mises stress at Node a 1 is reduced slightest and that at Node a 5 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 Figure6 b).The stress amplitudes of Nodes a 1 and a 5 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 a 1 or Node a 5 .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.Figure7shows 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 figure6, similar conclusions can be obtained.The stress amplitudes of Nodes b 1 and b 5 are larger than those of other three nodes, but smaller than those of Nodes a 1 and a 5 .Accordingly, under harmonic displacement, fatigue fracture is most likely to happen at the bottom of the first thread.

Figure 6 :
Figure 6: Von Mises Stress of the bottom of the first root of the thread.
Figure 7: Von Mises Stress of the edge of the contact area of the first root of the thread.

Figure 8 :
Figure 8: Variations of the maximum and minimum of the ratio of the friction stress to the shear stress with loading history.
Figure 11 shows the hysteresis loops for different bolt preloads.When a preload of 12.5 kN is used, friction forces between contacting surfaces are small.If rev -0.33 xx  mm, Contact I and Contact II are under the gross slip condition, and the shear load remains constant.With the increa-

Figure
Figure 10: Hysteresis loops for varying Coefficients of friction b  and