Dynamic Analysis of Bolted Cantilever Beam with Finite Element Method Considering Thread Relation

: There are complex nonlinear behaviors and mechanisms in the bolted joint interface. Thus, the bolted joint is crucial to the complex nonlinear dynamic response of the structure. However, in the traditional structural dynamic analysis, the screw connection is usually neglected, which makes it challenging to analyze and study the nonlinear dynamic behavior of bolted structures. Hence, based on the Timoshenko beam theory and ﬁnite element method, this paper introduces a model considering thread connection to analyze the dynamic response under di ﬀ erent excitation. Eventually, the results indicate that owing to the local nonlinearity of bolts, the whole bolted cantilever beam shows hardening-type characteristics. In addition, the frequency response curve also depicts the typical nonlinear phenomenon of instability and uncertainty, namely bifurcation, which preliminarily veriﬁes the correctness and accuracy of the bolted cantilever beam model.


Introduction
Bolted connection is extensively used in assembly structure, and the existence of interface leads to the discontinuity of local stiffness and damping.Particularly in the vibration environment, the relative slip of the interface may occur in the tangential direction, and the gap separation and impact collision may occur in the normal direction.The slip forms include micro slip, which only occurs in the local area of the interface and macro slip in the macro scale.The gap and collision forms include both micro slap and macro separation at the macro scale.Adhesion sliding and clearance collision are two typical nonlinear behaviors of bolted connections in vibration environment.Under the control of these behaviors, the dynamic response of the structure shows rich and complex nonlinear vibration phenomena, such as frequency doubling, frequency division, bifurcation, chaos and so on.
Tian [1] proposed to use materials with low elastic modulus and low Poisson's ratio to replace the description of joint surface.Thus, the nonlinear problem of structure was transformed into the nonlinear problem of material.Iranzad [2] exploited virtual elastoplastic materials to simulate the nonlinear behavior of the interface.The elastic modulus, Poisson's ratio, and shear modulus of the virtual material were obtained by taking the natural frequency of the modal obtained from the experimental test as the optimization result.Through harmonic response analysis, the macro slip force and stiffness of bolt joint surface were identified.Then, the virtual material model was defined and verified under different loads, which showed that the virtual material model could simulate the behavior of adhesion, micro slip, and macro slip of the joint surface.
According to Pohrt [3], based on the fractal contact theory, the boundary element method was used to analyze the normal stiffness and shrinkage resistance of two surfaces with arbitrary roughness and fractal dimension varying in a certain range.Jiang [4] used the fractal contact theory of rough surfaces to reveal the elastic-plastic deformation of contact surfaces, the size dependence of contact stiffness of micro contacts, and the source of contact stiffness in Hertz contact theory.Buczkowski [5] suggested a new fractal model based on the original assumption of contact area distribution based on the contact fractal model of single variable Weierstrass-Mandelbrot function, taking the roughness deformation and coupling effect into account and the results showed that the contact stiffness had a certain nonlinear evolution with the change of pressure.
Ahmadian [6] introduced Euler Bernoulli beam to establish the model of bolt connection structure.In the bolt connection, the stiffness was softened due to sliding.A nonlinear spring was used to replace the contact, and the stiffness and damping parameters of the joint surface are identified by harmonic response experiment, and then the elastic spring stiffness and damping were defined for the nonlinear dynamic modeling of the joint surface.Hanss [7] studies the bolt connection structure where the behavior of the bolt joint surface is replaced by spring and damping element.Through the test of force hammer and laser displacement sensor response transmission, the stiffness and damping parameters of the joint surface were identified based on fuzzy algorithm.Jalali [8] identified the parameters of the joint surface based on the force state mapping method, assuming that the response of acceleration and force in time domain was approximate to the single frequency excitation of the first natural frequency, so the modal stiffness and damping parameters were identified to represent the nonlinear behavior of the joint surface.Tsai [9] utilized the least square method to derive joint attributes directly from the frequency response functions measured in substructures and assembly structures, and the attributes extracted from one group of structures could also predict the frequency response functions of another group of structures with the same joint surface parameters, indicating that the method was effective.
Cescoto [10] advised an original method for numerical simulation of unilateral contact using finite element method.The main idea was that the displacement field and the contact stress field (pressure and friction shear) were discretized independently.Based on the variation principle, the theory was first proposed in the framework of infinitesimal deformation, and then extended to large inelastic strain.Lehnhoff [11,12] exploited contact finite element method to establish bolt connection model and calculated the displacement and contact pressure distribution of bolt joint surface under different bolt size, different thickness of connected parts and different material conditions, then gave stiffness curve and empirical formula of several bolt connection structures.Krishna [13] employed contact finite element algorithm to calculate the effect of joint in flange structural seal.Through experiments and finite element analysis of gasket loading and unloading characteristics, it was found that the contribution of the stress distribution to the sealing performance was greater than that of the flange rotation.Nassar [14] apply the finite element method to calculate the separation force of bolt connection structure under the action of separation force, which was distributed in the position of 3-5 times of bolt diameter, and other parts were in the state of no extrusion and no force transmission.Kim [15] studied four finite element modeling methods for bolted joint surfaces, which were solid element and contact element, coupling degree of freedom, beam element acting as bolt, and bolt free model with preload only applied between two plates.Modal test results showed that the error between the model established by solid element and contact element and experimental results were the least.Chen [16] calculated the influence of friction, clearance, bolt elasticity, stacking sequence and preload on the pressure distribution of multi-layer structure based on contact finite element algorithm, and carried out experimental verification.McCarthy [17][18][19] established a three-dimensional contact finite element model and experimental study, which proved that the pressure distribution of bolt plate connection structure was almost not affected by the gap between bolt and hole.With the increase of the clearance, the rotation angle of the bolt increases, the contact area decreased, and the contact stiffness increased.Ju [20] established the contact model of the joint surface to study the nonlinear characteristics of the bolted structure.It showed that the nominal force obtained from the finite element analysis was almost linearly proportional to the number of bolts arranged in the connection.Kucharski [21] inaugurate d a contact model of a micro hemisphere and a flat plate through finite element method and explored the micro mechanism of the joint surface.
Iwan [22,23] raised an Iwan model in which numerous Jenkins elements (spring and critical force elements) were connected in parallel to simulate the hysteretic characteristics of structures.Segalman [24] and Li [25] created four and six parameter continuous Iwan model, respectively, and identified the undetermined parameters according to the hysteretic curve of restoring force and displacement, thus revealing the relationship between energy dissipation and friction and preload.Brake [26] computed the Iwan model with four parameters and five parameters and discussed how to use other distribution functions to describe the distribution of critical force in Iwan model.Based on the theory of continuous Iwan model and beam element, Song [27,28] skillfully transformed the solution of Iwan model parameters into the calculation of transverse and torsional displacement in Rayleigh damping by using neural network algorithm.The Rayleigh damping parameters of bolted structures were identified by three-layer neural network algorithm, and the dynamic response is verified, which shows the effectiveness of the method.Segalman [29] investigated a three legged structure with bolt connection.Jenkins element was used to model in transverse and torsional directions at the joint, and the degrees of freedom in other directions were simplified as linear springs.The comparison of the response of the structure under different frequencies with that of the simplified model proved the effectiveness of the simplified model.Miller [30] scrutinized Iwan model to simulate the tangential hysteretic characteristics of the joint surface and used the model to simplify the model.The validity of the model was verified by response transfer experiment and response analysis.
Liu [31] created a three-dimensional finite element model in the ABAQUS to simulate the bolted joint under axial excitation, and studied the frictional stress, slip amplitude, and frictional work per unit area along two specified paths on the first thread.It was found that the finite elementresults agree with the experimental observations superbly.The Iwan model, however, was not extensively adopted owing to the high computational expense in a numerical simulation.Romagnoli [32] formulated a new experimentally validated expression for the bearing friction torque component in threaded joints with spherical and conical contact geometry, by extending two-dimensional Motosh model for various scenarios of contact pressure.In order to improve the understanding and the development of pertinent models, Daouk [33] designed a dynamic test bed, based on a bolted structure, and performed modal testing.The results of the experimental campaign showed the variation of the dissipation in a bolted joint and its apparent stiffness was a function of joint conditions.Adel [34] investigated the FE bolted joints model updating of a hybrid metal/composite structure and proposed a new approach called "doubly connective layer" for hybrid interfaces.The updated model compared to the initial model simulates experimental results more appropriately.Saito [35] discussed the dynamic characteristics of plastic plates with bolted joints and investigated the effects of tightening torque on the modal properties of the plates.By comparing the experimental and the numerical results, a modeling strategy for the boundary conditions between the plates was introduced and its validity was verified.
On the basis of literature, most researchers mainly focus on the local nonlinearity caused by the nonlinear spring derived from the nonlinear bolt connection model.The stable and unstable bifurcation phenomena of nonlinear systems with local and nonlocal nonlinearities are studied.In addition, in the previous work, most researchers used the Euler Bernoulli model assumption to simplify the lateral vibration of bolted structures.This theory ignores the influence of shear deformation and moment of inertia on the dynamic response of bolted structures.Therefore, based on Timoshenko beam theory and finite element method, this paper introduces the model of screw connection and analyzes the dynamic response under different excitation.The results show that due to the local nonlinearity of the bolt, the whole bolt cantilever beam presents the characteristics of hardening type.In addition, the frequency response curve also displays the typical nonlinear phenomenon of instability and uncertainty, namely bifurcation, which preliminarily verifies the correctness and accuracy of the bolt cantilever beam model.

Materials and Methods
Figure 1 depicts the considered bolted cantilever beam that consists of one cantilever beam in the left and a free lap beam in the right connected by a bolt.It is assumed that one single sine excitation force is applied in the right end of the structure, with being the excitation amplitude and being the excitation frequency.According to Figure 1, the x, y, and z represent the axial, lateral and transverse direction, respectively.Other parameters of the bolted cantilever beam, s, l, h, and L denote the distance from the central axis of the bolt to the wall, the length and thickness of a single beam, the height of beam, and the length of the whole structure of the bolt connected beam, respectively.Commonly, scholars only simplify the bolt to the transverse and normal contact problems, and do not consider the physical characteristics of the bolt.Bolt connection is different from mortise joint and riveting joint.It has strict connection mode of thread pair.Therefore, a finite element model considering thread characteristics is established in this paper and dynamic response analysis is conducted under different external excitation.Screw pair is composed of screw and nut.It is a kind of spatial motion pair, and its contact surface is spiral surface.When the screw and nut threads are subjected to axial load, when the screw or nut is turned, the relative sliding will occur between the two threads, and there is friction between the screw surfaces.When studying the friction in a screw pair, it is usually assumed that the axial load v G between the screw and the nut acts on the helix with pitch diameter.Since the helix can be developed into an oblique straight line on the plane, the effect of the force in the screw pair is the same as that between the slider and the inclined plane, as shown in Figure 3.In this way, the friction Figure 2 is a sectional view of the interaction between the screw and the nut.According to Figure 2, G v represents the axial force on the screw, β v represents the half angle of thread profile and θ v represents the half angle of groove.Obviously, the relationship between half angle of thread profile and half angle of groove is as follows, θ v = 90 • − β v .
Appl.Sci.2020, 10, x FOR PEER REVIEW 4 of 18 denote the distance from the central axis of the bolt to the wall, the length and thickness of a single beam, the height of beam, and the length of the whole structure of the bolt connected beam, respectively.Screw pair is composed of screw and nut.It is a kind of spatial motion pair, and its contact surface is spiral surface.When the screw and nut threads are subjected to axial load, when the screw or nut is turned, the relative sliding will occur between the two threads, and there is friction between the screw surfaces.When studying the friction in a screw pair, it is usually assumed that the axial load v G between the screw and the nut acts on the helix with pitch diameter.Since the helix can be developed into an oblique straight line on the plane, the effect of the force in the screw pair is the same as that between the slider and the inclined plane, as shown in Figure 3.In this way, the friction problem in the spatial screw pair can be transformed into the friction problem in the grooved moving Screw pair is composed of screw and nut.It is a kind of spatial motion pair, and its contact surface is spiral surface.When the screw and nut threads are subjected to axial load, when the screw or nut is turned, the relative sliding will occur between the two threads, and there is friction between the screw surfaces.When studying the friction in a screw pair, it is usually assumed that the axial load G v between the screw and the nut acts on the helix with pitch diameter.Since the helix can be developed into an oblique straight line on the plane, the effect of the force in the screw pair is the same as that between the slider and the inclined plane, as shown in Figure 3.In this way, the friction problem in the spatial screw pair can be transformed into the friction problem in the grooved moving pair.
where h P is the thread lead, n is the numbers of thread heads, v d is pitch diameter and P is the pitch of threaded bolt.
When tightening the nut, it is equivalent to the sliding block rising along the inclined plane.The peripheral force applied on the pitch diameter of the nut can be calculated as tan( ) The moment of the circumferential force to the thread axis is equal to the tightening torque required when tightening the nut; thus, tan( ) 22 Similarly, as is shown in Figure 4, when loosening the nut, it is equivalent to the sliding block sliding along the inclined plane at the same sliding speed, and the circumferential force applied on the pitch diameter of the nut can be calculated as follows: tan( ) Figure 3 shows a triangular screw pair with axial force and G v acting on the nut (including the weight of the nut itself).F R is the reaction force of screw on nut, and F is the circumferential force acting on thread pitch circle.Since the nut can be simplified as a slider, and the screw thread can develop into an inclined groove surface with an inclination angle λ v of along the middle diameter, the relationship between the nut and the screw can be simplified as that of the slider sliding along the groove surface.The inclination angle of the inclined plane is the thread rising angle on the pitch diameter of the thread; thus, where P h is the thread lead, n is the numbers of thread heads, d v is pitch diameter and P is the pitch of threaded bolt.When tightening the nut, it is equivalent to the sliding block rising along the inclined plane.The peripheral force applied on the pitch diameter of the nut can be calculated as The moment of the circumferential force to the thread axis is equal to the tightening torque required when tightening the nut; thus, Similarly, as is shown in Figure 4, when loosening the nut, it is equivalent to the sliding block sliding along the inclined plane at the same sliding speed, and the circumferential force applied on the pitch diameter of the nut can be calculated as follows: Appl.Sci.2020, 10, x FOR PEER REVIEW 6 of 18 The moment of the circumferential force M  to the thread axis is the tightening torque required to loosen the nut; thus, tan( ) 22 When vv   , M  is the resistance torque to prevent the nut from accelerating loosening while vv   , M  is the driving torque required to make the nut loose at the same speed.According to Figures 1-4, in order to model the complex interface of bolted joint, a novel finite element model of the bolted cantilever beam can be obtained in  The total kinetic energy of the beam system incorporates the kinetic energies of bolted cantilever beam and free beam, and they are shown in the below: where  , A , and I are density, cross section, and area moment of inertia of both of beams, respectively.The moment of the circumferential force M to the thread axis is the tightening torque required to loosen the nut; thus, When λ v > ϕ v , M is the resistance torque to prevent the nut from accelerating loosening while λ v < ϕ v , M is the driving torque required to make the nut loose at the same speed.
According to Figures 1-4, in order to model the complex interface of bolted joint, a novel finite element model of the bolted cantilever beam can be obtained in Figure 5.The bolted cantilever consists of a linear translational viscous damping c l , a nonlinear translational spring k ln , a linear torsional viscous damping c θ and a nonlinear torsional spring k θn .Moreover, uniform mass distribution is assumed to dodge the intricacy in the theoretical process in this paper.In addition, u(x, t), w(x, t), and φ(x, t) are employed to denote the axial, transverse, and rotation displacements as shown in Figure 6.The moment of the circumferential force M  to the thread axis is the tightening torque required to loosen the nut; thus, tan( ) 22 When vv   , M  is the resistance torque to prevent the nut from accelerating loosening while vv   , M  is the driving torque required to make the nut loose at the same speed.According to Figures 1-4, in order to model the complex interface of bolted joint, a novel finite element model of the bolted cantilever beam can be obtained in Figure 5  The total kinetic energy of the beam system incorporates the kinetic energies of bolted cantilever beam and free beam, and they are shown in the below: where  , A , and I are density, cross section, and area moment of inertia of both of beams, respectively.The total kinetic energy of the beam system incorporates the kinetic energies of bolted cantilever beam and free beam, and they are shown in the below: where ρ, A, and I are density, cross section, and area moment of inertia of both of beams, respectively.On the other hand, the total potential energy of the bolted cantilever beam system, which comprise potential energy of both of the beams, nonlinear spring, and nonlinear torsional spring as shown below, where U b1 and U b2 are potential energies of the cantilever beam and the lap beam, respectively.U nl and U θn are potential energies of the nonlinear translational spring and the nonlinear torsional spring, respectively.
According to Timoshenko beam theory, any axial displacement which can be caused by shear deformation is neglected in present exploration; thus, where w is the total transverse displacement of the bolted cantilever beam.w s and w b are shear and bending displacements of the bolted cantilever beam, respectively.β denotes the shear angle the bolted cantilever beam, and v is the lateral displacement the bolted cantilever beam.Thus, the strains corresponding to the displacement fields can be obtained as follows, For the sake of completeness, the effects of the axial displacement in the equation of motion of the bolted cantilever beam structure are taken into account.So, the component of stress corresponding to the strains of Equation ( 13) are given by (13) where E, G, and k are Young's modulus, shear modulus and shear factor for the rectangle cross section, correspondingly.
According to theory of solid mechanics, the strain energy of beam U b can be printed as By substituting Equations ( 12) and ( 13) into Equation ( 14), the strain energy of bolted cantilever beam can be obtained as follows: where I and A are moment of inertia and cross section area of the both of the beams.The potential energies of nonlinear spring and nonlinear torsional spring can be printed as follows: Appl.Sci.2020, 10, 7036 8 of 16 The variation of work done by external forces, which is generally comprise of the sine excitation force and viscose damping force from the surrounding, and these external forces are expressed as where c nl and c θn are damping coefficient of the bolted joint in translation and rotation directions of the cantilever beam system, respectively, and c is the viscos damping coefficient of the surroundings.
Then, the variation principle is applied to derive the nonlinear equation of the nonlinear dynamic for bolted cantilever beam.
where δ indicates variation operator.Substituting Equations ( 8)-( 10), (18), and ( 19) into Equation ( 20), nonlinear differential equations of the dynamic system of the bolted cantilever beam can be gained as follows: Although the dimensionality index is more sensitive to the signal characteristics, it will also change with the change of working conditions, and is easily affected by environmental interference, which has the defect of unstable performance.In contrast, the dimensionless index can eliminate the influence of these disturbance factors.So, the following dimensionless parameters in order to find the dimensionless nonlinear equations of the dynamic system of the bolted cantilever beam are presented: Substituting parameters introduced in Equation ( 27) into Equations ( 21)-( 26), and dropping asterisk notation for briefness.Thus, nonlinear differential equation of the dynamic structure can be translated into dimensionless partial differential equations.The results are as follows: Appl.Sci.2020, 10, 7036 9 of 16 To obtain the reduced order model of the partial equations of the dynamic system, the Galerkin method is introduced to discretize the whole structure, and the approximate series expansion is as follows: where ϕ k (x) and ψ k (x) are the k th Eigen function for the transverse and bending vibrations of the whole bolted cantilever beam, respectively.In order to obtain a couple of nonlinear ordinary differential equations of the bolted cantilever beam, Equation ( 34) is introduced into Equations ( 28)- (33), and the resultant equations are multiplied by the Engen function and integrated with x from 0 to 1.The nonlinear ordinary differential equations are shown as follows: .. .. .. .. ..
where q ij , p ij , and r ij are representing transverse displacement, bending displacement, and axial displacement in the generalized coordinate, respectively.Where N i means coefficient relating to the shape function of every beam.
For the convenience of reading, all of the important symbols have been given in the Nomenclature.
Figure 6a,b describe that the peaks of the frequency response curve are inclined to the right.This asymmetry reflects the existence of a certain degree of nonlinearity in transverse and bending vibration.Specifically, near the right side of the resonance frequency, there is a frequency corresponding to two amplitude of vibration, which is to say, there is a certain degree of uncertainty.Furthermore, no matter the transverse vibration or the bending vibration, there will be typical bifurcation phenomenon.Of course, if the bifurcation phenomenon is further developed, it may cause strong uncertainty chaos.Figure 7a,b show that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p initially slowly rise to their respective resonance peaks, and then rapidly drop to around 0. All of them will display a certain scale of bifurcation.
Moreover, with the increase of  , the resonance frequencies in both transverse direction and rotation direction increase in turn, and the formant decrease in turn.It may be that the natural frequency will increase as soon as it is too large.Under the action of small external force, there is not

Results and Discussion
Physical parameters of bolted cantilever beam are listed in Table 1, and the numerical results are presented in Figures 6-11.It can be seen from Figure 6a that in the first order generalized coordinate system, with the dimensionless frequency Ω increasing from 1 to 2, the generalized transverse coordinate value q slowly increases from 8.2 × 10 −7 , and the maximum value is 1.402 × 10 −5 , and the corresponding dimensionless resonance frequency is 1.642.Then, the value of the generalized horizontal coordinate decreased rapidly to 4.827 × 10 −7 .
Appl.Sci.2020, 10, x FOR PEER REVIEW 12 of 18 enough time for the bolt connected cantilever system to make response, so as to increase the amplitude.Figure 8a,b exhibit that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p primarily slowly rise to their respective resonance summits, and then rapidly drop to about 0. All of them will exhibit a certain level of bifurcation.In addition, with the gains of  , the resonance frequencies increase in sequence in transverse direction or rotation direction, and the formant decrease in turn.It may be that if the mass is too large, the natural frequency will be reduced.If the mass is too large, the inertia will increase.Under the action of small external force, the original vibration state of bolt connected cantilever beam system is not easy to change, so the amplitude becomes smaller.Figure 10a,b indicate that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p first gently upswing to their respective resonance peaks, and then speedily drop to near 0. All of them will reveal a certain amount of bifurcation.With the increase of ln k and nonlinear displacement stiffness, the natural frequency resonance frequencies of the dynamic bolted system increased in both transverse direction and rotation direction, and the maximum amplitude of the system is reduced.It may be that the natural frequency will increase as soon as it is too large.Under the action of external force, there is not enough time to increase the amplitude of the bolt connected cantilever beam system.Figure 10a,b indicate that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p first gently upswing to their respective resonance peaks, and then speedily drop to near 0. All of them will reveal a certain amount of bifurcation.With the increase of ln k and nonlinear displacement stiffness, the natural frequency resonance frequencies of the dynamic bolted system increased in both transverse direction and rotation direction, and the maximum amplitude of the system is reduced.It may be that the natural frequency will increase as soon as it is too large.Under the action of external force, there is not enough time to increase the amplitude of the bolt connected cantilever beam system.Figure 6b displays that in the first order generalized coordinate system, with the dimensionless frequency Ω increasing from 1 to 2, the generalized bending coordinate value p slowly enlarges from 2.216 × 10 −6 to a maximum value of 2.853 × 10 −5 , and the corresponding dimensionless resonance frequency is 1.660.Then, the value of the generalized horizontal coordinate decreased rapidly to 6.846 × 10 −7 .Figure 12 depicts that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases, the generalized transverse coordinate values q slowly upsurge to their respective resonance peaks, and then rapidly drop to near 0. All of them will expose a certain level of bifurcation and hard type characteristics.Compared with the Euler Bernoulli beam, the first order generalized coodinate of the frequency response curve of the Timoshenko beam in transverse direction has similar behavior, but amplitude and resonance frequency are lower than that of the Euler Bernoulli beam.The shear deformation is considered in the Timoshenko beam, which results in the internal energy dissipation due to the shear deformation, so the resonance peak of the Timoshenko beam is slightly lower than that of the Euler Bernoulli beam.Additionally, due to the influence of the moment of inertia, the stiffness of the Timoshenko beam becomes smaller and the resonance frequency becomes smaller.Figure 6a,b describe that the peaks of the frequency response curve are inclined to the right.This asymmetry reflects the existence of a certain degree of nonlinearity in transverse and bending vibration.Specifically, near the right side of the resonance frequency, there is a frequency corresponding to two amplitude of vibration, which is to say, there is a certain degree of uncertainty.Furthermore, no matter the transverse vibration or the bending vibration, there will be typical bifurcation phenomenon.Of course, if the bifurcation phenomenon is further developed, it may cause strong uncertainty chaos.
Figure 7a,b show that in the first-order generalized coordinate system, as the value of dimensionless frequency Ω increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p initially slowly rise to their respective resonance peaks, and then rapidly drop to around 0. All of them will display a certain scale of bifurcation.Moreover, with the increase of β, the resonance frequencies in both transverse direction and rotation direction increase in turn, and the formant decrease in turn.It may be that the natural frequency will increase as soon as it is too large.Under the action of small external force, there is not enough time for the bolt connected cantilever system to make response, so as to increase the amplitude.
Figure 8a,b exhibit that in the first-order generalized coordinate system, as the value of dimensionless frequency Ω increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p primarily slowly rise to their respective resonance summits, and then rapidly drop to about 0. All of them will exhibit a certain level of bifurcation.In addition, with the gains of η, the resonance frequencies increase in sequence in transverse direction or rotation direction, and the formant decrease in turn.It may be that if the mass is too large, the natural frequency will be reduced.If the mass is too large, the inertia will increase.Under the action of small external force, the original vibration state of bolt connected cantilever beam system is not easy to change, so the amplitude becomes smaller.
Figure 9a,b reveal that in the first-order generalized coordinate system, as the value of dimensionless frequency Ω increases from 1 to 2, the generalized transverse coordinate value q and the generalized rotation coordinate values p first slowly rise to their respective resonance crest, and then swiftly fall to around 0. All of them present an extent degree of bifurcation.With the increase of f 0 , the resonance frequency basically in transverse direction or rotation direction do not vary, which indicates that the vibration amplitude has little influence on the resonance frequency of the bolt connected beam.
Figure 10a,b indicate that in the first-order generalized coordinate system, as the value of dimensionless frequency Ω increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p first gently upswing to their respective resonance peaks, and then speedily drop to near 0. All of them will reveal a certain amount of bifurcation.With the increase of k ln and nonlinear displacement stiffness, the natural frequency and resonance frequencies of the dynamic bolted system increased in both transverse direction and rotation direction, and the maximum amplitude of the system is reduced.It may be that the natural frequency will increase as soon as it is too large.Under the action of external force, there is not enough time to increase the amplitude of the bolt connected cantilever beam system.
Figure 11a,b depict that in the first-order generalized coordinate system, as the value of dimensionless frequency Ω increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p formerly slowly upsurge to their respective resonance peaks, and then rapidly drop to near 0. All of them will expose a certain level of bifurcation.With the increase of k θn , the natural frequencies and resonance frequencies of the system change slightly in both transverse direction and rotation direction, and the vibration amplitudes decrease obviously.It may be that the torsional stiffness does not affect the resonance frequency of the system but affects the natural mode of the vibration system.
Figure 12 depicts that in the first-order generalized coordinate system, as the value of dimensionless frequency Ω increases, the generalized transverse coordinate values q slowly upsurge to their respective resonance peaks, and then rapidly drop to near 0. All of them will expose a certain level of bifurcation and hard type characteristics.Compared with the Euler Bernoulli beam, the first order generalized coodinate of the frequency response curve of the Timoshenko beam in transverse direction has similar behavior, but amplitude and resonance frequency are lower than that of the Euler Bernoulli beam.The shear deformation is considered in the Timoshenko beam, which results in the internal energy dissipation due to the shear deformation, so the resonance peak of the Timoshenko beam is slightly lower than that of the Euler Bernoulli beam.Additionally, due to the influence of the moment of inertia, the stiffness of the Timoshenko beam becomes smaller and the resonance frequency becomes smaller.Figure 12 depicts that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases, the generalized transverse coordinate values q slowly upsurge to their respective resonance peaks, and then rapidly drop to near 0. All of them will expose a certain level of bifurcation and hard type characteristics.Compared with the Euler Bernoulli beam, the first order generalized coodinate of the frequency response curve of the Timoshenko beam in transverse direction has similar behavior, but amplitude and resonance frequency are lower than that of the Euler Bernoulli beam.The shear deformation is considered in the Timoshenko beam, which results in the internal energy dissipation due to the shear deformation, so the resonance peak of the Timoshenko beam is slightly lower than that of the Euler Bernoulli beam.Additionally, due to the influence of the moment of inertia, the stiffness of the Timoshenko beam becomes smaller and the resonance frequency becomes smaller.

Conclusions
Based on Timoshenko beam theory and finite element principle, this paper introduces the model of screw connection and analyzes the dynamic response under distinctive excitation forces.By comparing with the Euler Bernoulli beam, the first order generalized coordinate of the frequency response curve of the Timoshenko beam in transverse direction has similar behavior, but amplitude and resonance frequency are lower than that of the Euler Bernoulli beam.The influences of other parameters such as β, η, f 0 , k ln , and k θn in the first-order generalized coordinate of the frequency response curve are discussed as well.Moreover, the corresponding results exhibit that owing to the local nonlinearity of the bolted joint, the whole bolt cantilever beam system presents the characteristics of hardening type.Moreover, the frequency response curves also display the typical nonlinear phenomenon of instability and uncertainty, namely bifurcation, which preliminarily verifies the correctness and accurateness of the bolted cantilever beam model.The viscos damping coefficient of the surroundings δ variation operator ϕ k (x) The k th Eigen function for the transverse vibrations of the whole bolted cantilever beam ψ k (x) The k th Eigen function for the bending vibrations of the whole bolted cantilever beam q ij The transverse displacement in the generalized coordinate p ij The bending displacement in the generalized coordinate r ij The axial displacement in the generalized coordinate N i The coefficient relating to the shape function of every beam k nl * The normalized nonlinear translational spring coefficient k θn * The normalized nonlinear torsional spring coefficient t * The normalized time Ω The normalized frequency Appl.Sci.2020, 10, x FOR PEER REVIEW 4 of 18denote the distance from the central axis of the bolt to the wall, the length and thickness of a single beam, the height of beam, and the length of the whole structure of the bolt connected beam, respectively.

Figure 2
is a sectional view of the interaction between the screw and the nut.According to Figure2, vG represents the axial force on the screw, v  represents the half angle of thread profile and v  represents the half angle of groove.Obviously, the relationship between half angle of thread profile and half angle of groove is as follows,

Figure 1 .
Figure 1.Structural diagram of bolted cantilever beam.Commonly, scholars only simplify the bolt to the transverse and normal contact problems, and do not consider the physical characteristics of the bolt.Bolt connection is different from mortise joint and riveting joint.It has strict connection mode of thread pair.Therefore, a finite element model considering thread characteristics is established in this paper and dynamic response analysis is conducted under different external excitation.Figure2is a sectional view of the interaction between the screw and the nut.According to Figure2, G v represents the axial force on the screw, β v represents the half angle of thread profile and θ v represents the half angle of groove.Obviously, the relationship between half angle of thread profile and half angle of groove is as follows, θ v = 90• − β v .

Figure 2
is a sectional view of the interaction between the screw and the nut.According to Figure2, vG represents the axial force on the screw, v  represents the half angle of thread profile and v  represents the half angle of groove.Obviously, the relationship between half angle of thread profile and half angle of groove is as follows,

Figure 3 .
Figure 3. Kinematic relationship between screw and nut during tightening.

Figure 3 .
Figure 3. Kinematic relationship between screw and nut during tightening.Assume the half angle of the thread being β v = 30 • ; thus, the groove semi-angle is θ v = 60 • , and the corresponding equivalent friction coefficient and the equivalent friction angle of the groove semiangle are shown below separately,

Figure 4 .
Figure 4. Kinematic relationship between screw and nut during loosening.

Figure 5 .
The bolted cantilever consists of a linear translational viscous damping l c , a nonlinear translational spring ln k , a linear torsional viscous damping c  and a nonlinear torsional spring n k  .Moreover, uniform mass distribution is assumed to dodge the intricacy in the theoretical process in this paper.In addition, ( , ) u x t , ( , ) w x t , and ( , ) xt  are employed to denote the axial, transverse, and rotation displacements as shown in Figure 6.

Figure 4 .
Figure 4. Kinematic relationship between screw and nut during loosening.

Figure 4 .
Figure 4. Kinematic relationship between screw and nut during loosening.
. The bolted cantilever consists of a linear translational viscous damping l c , a nonlinear translational spring ln k , a linear torsional viscous damping c  and a nonlinear torsional spring n k  .Moreover, uniform mass distribution is assumed to dodge the intricacy in the theoretical process in this paper.In addition, ( , ) u x t , ( , ) w x t , and ( , ) xt  are employed to denote the axial, transverse, and rotation displacements as shown in Figure 6.

Figure 6 .
Figure 6.1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 6 .
Figure 6.1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 7 .
Figure 7. Influence of  in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 7 .
Figure 7. Influence of β in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 8 .
Figure 8. Influence of  in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure
Figure9a,b reveal that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases from 1 to 2, the generalized transverse coordinate value q and the generalized rotation coordinate values p first slowly rise to their respective resonance crest, and then swiftly fall to around 0. All of them present an extent degree of bifurcation.With the increase of 0 f , the resonance frequency basically in transverse direction or rotation direction do not vary, which indicates that the vibration amplitude has little influence on the resonance frequency of the bolt connected beam.

Figure 8 . 18 Figure 9 .
Figure 8. Influence of η in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.Appl.Sci.2020, 10, x FOR PEER REVIEW 13 of 18

Figure 10 .
Figure 10.Influence of ln k in 1st order generalized coodinate of the frequency response curve (a) in

Figure 9 . 18 Figure 9 .
Figure 9. Influence of f 0 in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 10 .
Figure 10.Influence of ln k in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure
Figure11a,b depict that in the first-order generalized coordinate system, as the value of dimensionless frequency  increases from 1 to 2, the generalized transverse coordinate values q and the generalized rotation coordinate values p formerly slowly upsurge to their respective resonance peaks, and then rapidly drop to near 0. All of them will expose a certain level of bifurcation.With the increase of n k  , the natural frequencies and resonance frequencies of the system change slightly in both transverse direction and rotation direction, and the vibration amplitudes decrease obviously.It may be that the torsional stiffness does not affect the resonance frequency of the system but affects the natural mode of the vibration system.

Figure 10 .
Figure 10.Influence of k ln in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 12 .
Figure 12.Comparison between Timoshenko beam and Euler Bernoulli beam in 1st order generalized coodinate of the frequency response curve.
Based on Timoshenko beam theory and finite principle, this paper introduces the model of screw connection and analyzes the dynamic response under distinctive excitation forces.By comparing with the Euler Bernoulli beam, the first order generalized coordinate of the frequency response curve of the Timoshenko beam in transverse direction has similar behavior, but amplitude and resonance frequency are lower than that of the Euler Bernoulli beam.The influences of other parameters such as  ,  , 0 f , ln k , and n k  in the first-order generalized coordinate of the frequency response curve are discussed as well.Moreover, the corresponding results exhibit that

Figure 11 .
Figure 11.Influence of k θn in 1st order generalized coodinate of the frequency response curve (a) in transvese direction (b) in rotation direction.

Figure 12 .
Figure 12.Comparison between Timoshenko beam and Euler Bernoulli beam in 1st order generalized coodinate of the frequency response curve.
Based on Timoshenko beam theory and finite element principle, this paper introduces the model of screw connection and analyzes the dynamic response under distinctive excitation forces.By comparing with the Euler Bernoulli beam, the first order generalized coordinate of the frequency response curve of the Timoshenko beam in transverse direction has similar behavior, but amplitude and resonance frequency are lower than that of the Euler Bernoulli beam.The influences of other parameters such as  ,  , 0 f , ln k , and n k  in the first-order generalized coordinate of the frequency response curve are discussed as well.Moreover, the corresponding results exhibit that

Figure 12 .
Figure 12.Comparison between Timoshenko beam and Euler Bernoulli beam in 1st order generalized coodinate of the frequency response curve.
c nlNonlinear damping coefficient of the bolted joint in translation direction c θn Nonlinear damping coefficient of the bolted joint in rotation direction c

Table 1 .
Physical parameters of bolted cantilever beam.