Dynamic and Stability Analysis of Multibolt Plane Joints under Normal Forces

In this paper, a stiffness model of contact surfaces based on a modified three-dimensional fractal contact model is built, which is in accordance with the experiment results. Additionally, the static, dynamic, and stable behaviors of the bolt joint between the spindle box and the machine bed are analyzed. The mathematical relationship between fractal parameters of the surface topography and the stiffness of the system was established to accurately study its static behaviors. Asymmetric curves are observed from the load–deflection results and the nonlinear stiffness characteristic is also presented. It is shown that both the stress and the stiffness increase with the increase of the displacement near the static equilibrium position. Meanwhile, a simplified model without the consideration of roughness is compared with joint interfaces composed from milling, scraping, and grinding surfaces. Numerical calculation was employed to investigate effects of design parameters on the system under harmonic excitation when the processing method, excitation force, bolt pre-tightening force, topography parameters, and other structural parameters, i.e., nominal contact area, joint thickness and bolt number, are eventually regarded as the control parameters. The aim of the article is to analysis the influence of these parameters, including surface morphology, on nonlinear characteristics of the bolt interface with fractal contact surfaces andto provide some references to improve the characteristics.


Introduction
The contact interfaces between parts are known as mechanical joint surfaces. These interfaces play an important role in the dynamic behavior of mechanisms and mechanical devices. For example, contact stiffness and damping of such systems at machine tools contribute 60% of the overall stiffness and more than 90% of the overall damping [1,2]. Furthermore, many studies have focused on these joint surfaces. Liu et al. analyzed the dynamic characteristics of the joints in the aero-engine rotor system through FEM simulation and experiment [3]. A virtual gradient material model was proposed by Liao et al. to obtain the contact pressure distribution and dynamic characteristics of the bolted joint [4]. Daouket al. researched the dynamic behavior of the bolted joint under heavy loadings utilizing a bolted structure test bed [5]. Pan et al. built a critical two-degree-of-freedom modal coupling model and studied the effect of contact stiffness on the stability and nonlinearity on the basis of fractal surface topography [6]. Jalaliet al. provided a new linear joint model and investigated the influence of preload and surface roughness quality on the identified equivalent stiffness [7]. A new lump model applied in finite element analysis was given by Beaudoin et al. for the nonlinear dynamic response of bolted flanges [8]. With these interfaces withstanding contact pressures and frictions, the contact problems in engineering include geometry errors, microscopic roughness, lubrication, coating, conduction of heat, and electricity, which are nonlinear and complicated. Therefore, studies of the dynamic characteristics of mechanical joint surfacesare of great significance.
Research on modeling iscrucial for developing simulations to analyze the dynamic motions of contact interfaces. Conventionally, Hertzian contact theory is utilized to model interacting surfaces. As Greenwood and Williamson consider the fact that all engineering surfaces are not absolutely smooth, they describe the contact interfaces with extended Hertzian theory [9]. Rough surfaces are approximated by a series of spherical summits which have the same radius and Gaussian height distributions. The method was first based on statistical representation of the surface topography and is closer to the actual conditions than previous models, and has an enormous influence on contact theory. In the following, many investigations were focused on the extension of G-W model by Bhushan et al. [10], McCool [11], Nayak [12], and Whitehouse et al. [13]. However, Majumdar and Bhushan as well as Sayles and Thomas suggest that statistical roughness parameters depend on both instrumental resolution and sampling length [14,15]. It was found that the Weierstrass-Mandelbrot (W-M) function satisfies all these properties of a surface profile [14,16]. Therefore, fractal theory has been employed to model engineering surfaces by Majumdar and Bhushan (M-B) due to their statistic self-affinity, which is scale-independent [14,[17][18][19]. Wang and Komvopoulos introduced the domain extension factor for microcontact size distribution to modify the M-B model [20]. Another fractal-based model was proposed by using a multivariate W-M function [21,22]. According to the previous model, Shuyun Jiang et al. deduced the relationship between stiffness and the real area of contact mathematically [23]. Xiao et al. analyzed the static force-displacement characteristic of joint surfaces on the basis of a two-variable W-M function using the finite element method and further the dynamic behavior of the contact system [22,24].
The above mentioned studies mostly focused on the experiment test or FE method, which is often considered as a macroscopic system. In addition, the previous studies rarely referred to chaos analysis. In this work, the lumped-mass method is applied for the sake of simplicity of computation and the fractal contact model is employed to the perspective of the surface morphology. Based on a single-freedom dynamic model of the bolt joint between the spindle box and the machine bed, which is a typical structure in a machine tool, the influence of design parameters and the load on forced vibration responses has been discussed in order to clarify nonlinear characteristics of the bolt interface. The dynamic model applies to bolted plane joint parts with larger mass. A contact interface regardless of asperities is presented for the purpose of comparison. Moreover, the response analysis can help avoid resonance and chaos by adjusting the parameters in the design process. The outline of this article is as follows. Section 2 outlines a mathematical relationship of the stiffness per unit area and the pressure between joint engineering surfaces, which is verified by the experiment. In Section 3, a single-freedom dynamic model is formulated, where the normal stiffness and damping are treated as functions of fractal parameters. In Section 4, the load-deflection and the stiffness-deflection functions are investigated based onthe kinetic equations obtained. In Section 5, the fourth-order Runge-Kutta method is utilized to study effects of design parameters on the dynamic behaviors and stability of the bolt contact interface. Finally, conclusions are presented concisely in the Section 6.

Surface Fractal Contact Model and Verification
The contact interface is an important nonlinearfactor of the bolt joint system. A three-dimensional fractal load-separation contact model modified with the domain extension factor was introduced and the contact stiffness can be obtained byderivationof the curve.
To characterize features of the rough surface topography assuming isotropy, such as continuity, nondifferentiability, and self-affinity, a modified W-M function is given by [23]  where z(x) is the fractal profile, L is the sample length, G is the fractal roughness parameter, D (1 < D < 2) is the fractal dimension, γ (γ >1) is the scaling parameter controls the density of the frequencies in the surface profile, and ϕ 1n is a random phase. In consideration of surface flatness and frequency distribution density, γ = 1.5 is suitable [21]. The approximate root-mean-square height is obtained as [23] In previous studies, when two flexible bodies with nominally flat surfaces contact with each other, the system is usually replaced by a rigid smooth surface contact with an equivalent elastic rough surface with reduced elastic modulus to make calculation convenient, where E 1 , E 2 and ν 1 , ν 2 are elastic moduli and Poisson ratios of original surfaces, respectively [23,25]. The stress is derived from the pressure of contact asperities. To simplify the problem, it is normal and reasonable to assume that all asperities have spherical summits with different radiuses [17,21].
For the simple case of elastic-fully plastic material, the normal elastic force p e at an asperity can be written as [23] p e (a ) = 2 where a = 2a e is the truncated microcontact area belonging to a hypothetical section of the undeformed asperity at the contact position doubles the real microcontact area in elastic deformations and a = a p in fully plastic deformations. The critical truncated area a c is given by [23] (4) in which b = (πq/2) 2 and q = 0.454 + 0.41ν, where ν and H are the Poisson ratio and the hardness of the softer material, respectively. It is determined that asperities deform elastically when a > a c .
In addition, the distribution function of the truncated microcontact area is expressed as [20] n(a ) = where a l is the largest truncated microcontact area and ψ (ψ > 1) is the domain extension factor. The value of ψ can be obtained by dichotomy method depends on a transcendental equation ψ (2−D)/2 − For light pressure of the interface, asperities are far enough to ignore the effect of their interactions, which leads to a real contact area much smaller than the nominal contact area. A reasonable approximation of the total elastic force P e can be written as [21] P e a l = a l a c p e (a )n(a )da .
Then, substituting Equations (3) and (5) into Equation (6) yields The truncated contact area is related to n(a ) as follows [21]: A r a l = a l 0 a n(a )da (8) Introducing Equation (5) into Equation (8) gives If the rough surface obeys the Gaussian distribution, the truncated contact area also can be expressed as a function of the rms height, the separation x s between the mean planes of the joint interface and the nominal contact area A a , which is formed as [17] where erfc( ) is the complimentary error function and x s = x s0 − ∆x s , in which x s0 is the initial separation and ∆x s is the decrease of the distance between the two surfaces. Hence, the relationship of the separation x s and the force P e is obtained mathematically. The normal contact stiffness can be acquired by taking the derivative of P e with respect to x s .
In addition, Jiang et al. have proposed a stiffness model based on the fractal theory and carried out an experimental study utilized cast iron, which is the same as the material of the contact parts in this paper [23]. Figure 1 performs the comparison of numerical simulations of the present model, Jiang's model, and experiment results of Jiang. It shows the theory curve has good coherence to the experimental data. Here, the pressure of the interface Q and the stiffness per unit area K unit are obtained as Q a l = P e a l A a (11) where the stiffness of contact surfaces K s (∆x s ) = dP e (∆x s )/d∆x s .

Dynamic Model of the Bolt Joint Interface
In this section, a model of the spindle box in contact with the machine bed is built to analyze the dynamic characteristic of the fixed joint parts. For simplicity, contact objects are considered to be a series of asperities according to the equivalent surface on base B, two elastic smooth bases A (with a rigid bottom plane), B connected by bolts and two rigid bodies A, B, which is performed in Figure  2a. It has been assumed that pressures on bases are evenly distributed. The stiffness of serial bulks is

Dynamic Model of the Bolt Joint Interface
In this section, a model of the spindle box in contact with the machine bed is built to analyze the dynamic characteristic of the fixed joint parts. For simplicity, contact objects are considered to be a series of asperities according to the equivalent surface on base B, two elastic smooth bases A (with a rigid bottom plane), B connected by bolts and two rigid bodies A, B, which is performed in Figure 2a. It has been assumed that pressures on bases are evenly distributed. The stiffness of serial bulks is defined as K bases = K 1 + K 2 , in which K 1 = E 1 A a /l 1 and K 2 = E 2 A a /l 2 , where l 1 and l 2 are the thickness of bulks respectively. Deformation occurs at those asperities leading to a contact stiffness K s . The total stiffness of bolts K bolts = mk bolt is determined by the type and number m. As demonstrated in Figure 2a, K s and K bases are in series, and in parallel with K bolts . On the basis of the structure, a single-degree-of-freedom model of the system is established in Figure 2b. The parameter c is the damping coefficient of the system. The distance variation ∆x between the most remote two surfaces belong to bases and the variation ∆x s are caused by a harmonic excitation Fsin(ωt) in this model. Initial distances, x pre = x + ∆x and x spre = x s + ∆x s , corresponding to these two variables are both due to the mass of the spindle box M and the pretension force of bolts F pre = m f pre , where the subscript pre denotes the system is static equilibrium and preloaded by bolts. Figure 1. Dependence of the stiffness per unit area on the pressure: (a) milling/milling joint, (b) scraping/scraping joint, (c) grinding/grinding joint.

Dynamic Model of the Bolt Joint Interface
In this section, a model of the spindle box in contact with the machine bed is built to analyze the dynamic characteristic of the fixed joint parts. For simplicity, contact objects are considered to be a series of asperities according to the equivalent surface on base B, two elastic smooth bases A (with a rigid bottom plane), B connected by bolts and two rigid bodies A, B, which is performed in Figure  2a. It has been assumed that pressures on bases are evenly distributed. The stiffness of serial bulks is defined as = +   In order to obtain the force between the interfaceundervibration condition, aphenomenon is introducedthat even though the large plastic deformation occurs during the load process, the unloading process is completely elastic [26]. Therefore, the stressbetweencontact surfaces of a mechanical joint under repeated loading can beconsidered as an elastic force only, which is acquired by Equation (7). The kinetic equation is given by ∆x > ∆x bpre (13) in which x spre satisfied the function P e x spre = Mg + F pre , and ∆x bpre = f pre /k bolt is the elongation of preloaded bolts. Moreover, P e x spre − ∆x s − P e x spre = K bases (∆x − ∆x s ) is known according to Figure 2b to deduce the function ∆x = f (∆x s ) = P e x spre − ∆x s − P e x spre /K bases + ∆x s . Substituting this relationship into Equation (13) yields Introducing the overall stiffness at static balance K npre = K bases K spre / K bases + K spre + K bolts , where K spre = K s (0) depends on the value of x s0 = x spre , the dimensionless equation of the vertical vibration model can be obtained by defining τ = ωt, The equation of motion is transformed into the dimensionless form where By contrast, if the roughness of contact surfaces is neglected, the dynamic model changes into an ideal model and ∆x s can be expressed in the form (16) in which ∆x apre = − Mg+F pre K bases . Equation (14) can be transformed into As ∆x s = 0 in a small neighborhood of ∆x= 0, the contact stiffness at the static equilibrium position K spre is infinite, and then K npre = K bases + K bolts is achieved. Introducing the parameter G = Mg/K npre , the dimensionless equation is as follows The force between bases is obtained as

Force-Displacement Characteristic of the Bolt Joint Interface
To analyze the dynamic performance of the system, the stiffness should be obtained. On the basis of Equation (14) in Section 3, kinetic equation without inertia force and damping force illustrates the relationship between the normal external load F n and ∆x under slow load rate, which is given by Appl. Sci. 2019, 9, 5521 7 of 16 Furthermore, with the ignorance of unevenness, Equation (17) can be simplified as Taking the derivative of F n with respect to ∆x, the total stiffness K n is written as When roughnesses are overlooked, the expression of K n is transformed into The influence of different manufacturing methods on the static behavior of the system is shown in Figure 3. As demonstrated in Figure 3a, the slope of curves increases except for that of ideal surfaces, indicating a hardening nonlinear behavior when ∆x > 0 and a softening nonlinear behavior when ∆x < 0. This phenomenon is more apparently for scraping/scraping joint which has a smoother equivalent surface in −0.5 × 10 −5 m < ∆x < 0.5 × 10 −5 m. With the decrease of the displacement in the negative direction, these curves are closer to each other. It is caused by the interface separation of the contact pairs. Figure 3b shows that the stiffness without the consideration of asperities is piecewise-linear over the whole displacement range and greater than others when ∆x > 0. All the curves intersect each other, and the crossing points are observed in −0.5 × 10 −5 m < ∆x < 0. The parameters of bolt joint surfaces are shown in Tables 1 and 2 [23].       Values of topography parameters are the experiment data in the study by Jiang et al., a subscript eq denotes that D eq and G eq are equivalent fractal parameters.

Results
This section is concerned with effects of the dimensionless parameter, excitation force F, and design parameters, i.e., processing method, equivalent fractal dimension D eq , equivalent fractal roughness G eq , nominal contact area A a , thickness of the joint l = l 1 + l 2 , bolt number m and bolt pre-tightening f pre of the system on the dynamic behaviors and stability in excitation frequency ω utilizing the fractal parameters of the milling/milling joint in Table 2. The fourth-order Runge-Kutta method was employed to calculate the numerical solutions. Equations (7), (9), (10), and (13) indicate that nine parameters are independent for the present system, which are M, c, K bolts , D eq , G eq , A a , F, ω, and ∆x bpre . Thus, ω and F are converted into proper dimensionless indexes ω and F. Parameters c, K bolts , and ∆x bpre have exhibited some correlations existed with ζ, M, D eq , G eq , A a , l, m, and f pre , respectively, in previous analysis, where ζ and M are assumed to remain constant values. All the parameters of the joint interface are still provided by Table 1.

Effect of Manufacturing Methodon Nonlinear Response
In general, the surface appearance changes the dynamic behavior of the joint interface. The distinction among milling, scraping, grinding, and smooth surfaces was studied. For the amplitude of excitation, force F is fixed as 2000 N, and results of numerical simulations are demonstrated in Figure 4, whichexhibits softening type nonlinearity, and there appears to be a jumping phenomenon except for smooth surfaces. With the roughness R a of the machined surface increases, the primary resonance frequency and the maximum amplitude of the system increase. It is also found the system has a tendency of occurrence of superharmonic resonance near ω = 0.5.

Effect of Excitation Force on Nonlinear Response
In this subsection, the gradual change between period motion and chaotic motion in different dimensionless excitation forces F are shown in Figure 5. It is observed that a 2T-period response occurs from frequency ratio ω = 1.96 to ω = 2.03 for

Effect of Excitation Force on Nonlinear Response
In this subsection, the gradual change between period motion and chaotic motion in different dimensionless excitation forces F are shown in Figure 5. It is observed that a 2T-period response occurs from frequency ratio ω = 1.96 to ω = 2.03 for F = 2 × 10 −6 m. With the increase of F, another region of 2T-period motion near ω = 0.63 and region of 3T-period motion near ω = 1.4 appear in succession when F = 4 × 10 −6 m and F = 6 × 10 −6 m, respectively. As F is further increased, chaotic behavior is visible gradually at several points before ω = 0.6. In the meantime, 2T-and 3T-period motion regions and chaos regions are enlarged, and the midpoints of these regions move towards the negative direction.

Effect of Excitation Force on Nonlinear Response
In this subsection, the gradual change between period motion and chaotic motion in different dimensionless excitation forces F are shown in Figure 5. It is observed that a 2T-period response occurs from frequency ratio ω = 1.96 to ω = 2.03 for ,respectively. As F is further increased, chaotic behavior is visible gradually at several points before ω = 0.6 . In the meantime, 2T-and 3T-period motion regions and chaos regions are enlarged, and the midpoints of these regions move towards the negative direction. For − = × 6 7 10 m F , a 3-D waterfall chart with the frequency ratio as control parameter is illustrated in Figure 6. With the increase of ω , there appears a jumping phenomenon in the excitation frequency component at ω = 0.61 . It can also be seen that super-and subharmonic vibrations occur in different regions. The system generates a subharmonic, which leads to the bifurcation phenomenon. When the frequency spectrum is continuous near ω = 0.34 and ω = 0.54 , the system performs chaotic motion, which is according with data in Figure 5. For F = 7 × 10 −6 m, a 3-D waterfall chart with the frequency ratio as control parameter is illustrated in Figure 6. With the increase of ω, there appears a jumping phenomenon in the excitation frequency component at ω = 0.61. It can also be seen that super-and subharmonic vibrations occur in different regions. The system generates a subharmonic, which leads to the bifurcation phenomenon. When the frequency spectrum is continuous near ω = 0.34 and ω = 0.54, the system performs chaotic motion, which is according with data in Figure 5.  Figure 7, the joint interface shows a 1T-period motion at ω = 0.33 . However, as ω is increased to 0.34, the 1T-period motion is replaced by chaotic motion, which is seen in Figure 8.
When ω reaches 0.35, the system transits from chaotic motion to 2T-period motion is shown in Figure 9. Since the frequency spectrum is continuous and a number of discrete points appear in the Poincare section, chaotic behavior is obviously visible. At ω = 0.4 , the motion transits from 2T-period motion back to 1T-period motion is exhibited in Figure 10. From Figures 7-10, it can be observed that superharmonics exist in each frequency spectrum. To explain more clearly, Figures 7-10 demonstrate the time histories, frequency spectrums, phase plane diagrams, and Poincare sections under four different excitation frequencies for F = 7 × 10 −6 m. In Figure 7, the joint interface shows a 1T-period motion at ω = 0.33. However, as ω is increased to 0.34, the 1T-period motion is replaced by chaotic motion, which is seen in Figure 8. When ω reaches 0.35, the system transits from chaotic motion to 2T-period motion is shown in Figure 9. Since the frequency spectrum is continuous and a number of discrete points appear in the Poincare section, chaotic behavior is obviously visible. At ω = 0.4, the motion transits from 2T-period motion back to 1T-period motion is exhibited in Figure 10. From Figures 7-10, it can be observed that superharmonics exist in each frequency spectrum. To explain more clearly, Figures 7-10 demonstrate the time histories, frequency spectrums, phase plane diagrams, and Poincare sections under four different excitation frequencies for − = × 6 7 10 m F . In Figure 7, the joint interface shows a 1T-period motion at ω = 0.33 . However, as ω is increased to 0.34, the 1T-period motion is replaced by chaotic motion, which is seen in Figure 8.
When ω reaches 0.35, the system transits from chaotic motion to 2T-period motion is shown in Figure 9. Since the frequency spectrum is continuous and a number of discrete points appear in the Poincare section, chaotic behavior is obviously visible. At ω = 0.4 , the motion transits from 2T-period motion back to 1T-period motion is exhibited in Figure 10. From Figures 7-10, it can be observed that superharmonics exist in each frequency spectrum.

Effect of Fractal Parameters of Surface Topography
According to Equation (1), the surface topography is determined from the fractal dimension D and the fractal roughness G. Figure 11 shows that the stiffness with respect to the displacement in

Effect of Fractal Parameters of Surface Topography
According to Equation (1), the surface topography is determined from the fractal dimension D and the fractal roughness G. Figure 11 shows that the stiffness with respect to the displacement in

Effect of Fractal Parameters of Surface Topography
According to Equation (1), the surface topography is determined from the fractal dimension D and the fractal roughness G. Figure 11 shows that the stiffness with respect to the displacement in different equivalent fractal parameters. All the curves are monotonically increasing and closer to each other in the negative direction. As the displacement is increased from −3 × 10 −5 m to 1 × 10 −5 m, the curves intersect each other when −0.5 × 10 −5 m < ∆x < 0 and the stiffness increases with a smoother surface topography when ∆x > 0, i.e., D eq increases or G eq decreases, are illustrated in Figure 11a,b, respectively. different equivalent fractal parameters. All the curves are monotonically increasing and closer to each other in the negative direction. As the displacement is increased from −3 × 10 −5 m to 1 × 10 −5 m, the curves intersect each other when −0.5 × 10 −5 m < Δx < 0 and the stiffness increases with a smoother surface topography when Δx > 0, i.e., Deqincreases or Geq decreases, are illustrated in Figure 11a,b, respectively. The effect of Deq on nonlinear response is compared with the one of Geq. Figure 12 depicts that with the increase of Deq or decrease of Geq, three different regions appear in succession on the bifurcation diagrams: 1T-period motion, 2T-period motion, and 1T-period motion. As noted above, critical values Deq = 1.18 and Geq = 10 −11 m can be found for − = × 6 0.9 10 m F and − = × 6 1.1 10 m F respectively. Doubling bifurcations both occur just before ω = 2 and midpoints of 2T-period motion regions both move to the negative direction in Figure 12a,b. When the excitation frequency is out of the region of the 2T-period response, the change in displacement is slight and the dynamic behavior becomes relatively stable.

Effect of Other Structural Parameters on Nonlinear Response
The apparent contact area, thickness, and number of bolts are key parameters in the design of the bolt joint. Influences of the three parameters on the vibration of the system were studied. As ω is increased in the neighborhood of 2, a transition from 1T-period motion to 2T-period motion, and finally back to 1T-period motion again can be observed under certain conditions, i.e.,  The effect of D eq on nonlinear response is compared with the one of G eq . Figure 12 depicts that with the increase of D eq or decrease of G eq , three different regions appear in succession on the bifurcation diagrams: 1T-period motion, 2T-period motion, and 1T-period motion. As noted above, critical values D eq = 1.18 and G eq = 10 −11 m can be found for F = 0.9 × 10 −6 m and F = 1.1 × 10 −6 m respectively. Doubling bifurcations both occur just before ω = 2 and midpoints of 2T-period motion regions both move to the negative direction in Figure 12a,b. When the excitation frequency is out of the region of the 2T-period response, the change in displacement is slight and the dynamic behavior becomes relatively stable. different equivalent fractal parameters. All the curves are monotonically increasing and closer to each other in the negative direction. As the displacement is increased from −3 × 10 −5 m to 1 × 10 −5 m, the curves intersect each other when −0.5 × 10 −5 m < Δx < 0 and the stiffness increases with a smoother surface topography when Δx > 0, i.e., Deqincreases or Geq decreases, are illustrated in Figure 11a,b, respectively.  Figure 12a,b. When the excitation frequency is out of the region of the 2T-period response, the change in displacement is slight and the dynamic behavior becomes relatively stable.

Effect of Other Structural Parameters on Nonlinear Response
The apparent contact area, thickness, and number of bolts are key parameters in the design of the bolt joint. Influences of the three parameters on the vibration of the system were studied. As ω is increased in the neighborhood of 2, a transition from 1T-period motion to 2T-period motion, and finally back to 1T-period motion again can be observed under certain conditions, i.e.,

Effect of Other Structural Parameters on Nonlinear Response
The apparent contact area, thickness, and number of bolts are key parameters in the design of the bolt joint. Influences of the three parameters on the vibration of the system were studied. As ω is increased in the neighborhood of 2, a transition from 1T-period motion to 2T-period motion, and finally back to 1T-period motion again can be observed under certain conditions, i.e., A a > 0.03m 2 for F = 10 −6 m, l > 0.1m for F = 0.59 × 10 −6 m and m < 18 for F = 0.8 × 10 −6 m, are illustrated in Figure 13a-c, respectively. With the increase of A a or decrease of m, the 2T-period motion regions are enlarged. However, the width of the region changes monotonously with the thickness l, which is presumed to be a convex function. That is, the bolted joints tend to be stable if A a decreases, m increases or l increases after a certain value. Moreover, midpoints of 2T-period motion regions are not sensitive to the apparent contact area, thickness, and number of bolts. enlarged. However, the width of the region changes monotonously with the thickness l, which is presumed to be a convex function. That is, the bolted joints tend to be stable if Aa decreases, m increases or l increases after a certain value. Moreover, midpoints of 2T-period motion regions are not sensitive to the apparent contact area, thickness, and number of bolts.

Effect of Bolt Pre-Tightening on Nonlinear Response
The pretension force of bolts performs a great influence on the mechanical capacity of the system. In order to analyze the effect of the bolt pre-tightening on dynamic behaviors, bifurcation diagrams are presented for fpre = 3, 1, 0.1, and 0 kN. As demonstrated in Figure 14a, there appears a 2T-period response from frequency ratio ω = 1.99 to ω = 2.01. When fpre is decreased to 1 kN, the 2T-period motion region becomes narrow can be seen in Figure 14b. Meanwhile, another region of 2T-period response and a region of 3T-period response occur near ω = 0.65 and ω = 0.71 , respectively. Figure 14c depicts that as fpre continues to increase, more bifurcation points separate the region of excitation frequency into a number of stages. In some stages, 4T-, 5T-, 7T-, 8T-, 15T-, 30T-, and 40T-period responses are observed, which shows the appearance of subharmonics. With the further decrease of fpre, chaotic behavior at ω = 0.23 is clearly visible in Figure 14d.

Effect of Bolt Pre-Tightening on Nonlinear Response
The pretension force of bolts performs a great influence on the mechanical capacity of the system. In order to analyze the effect of the bolt pre-tightening on dynamic behaviors, bifurcation diagrams are presented for f pre = 3, 1, 0.1, and 0 kN. As demonstrated in Figure 14a, there appears a 2T-period response from frequency ratio ω = 1.99 to ω = 2.01. When f pre is decreased to 1 kN, the 2T-period motion region becomes narrow can be seen in Figure 14b. Meanwhile, another region of 2T-period response and a region of 3T-period response occur near ω = 0.65 and ω = 0.71, respectively. Figure 14c depicts that as f pre continues to increase, more bifurcation points separate the region of excitation frequency into a number of stages. In some stages, 4T-, 5T-, 7T-, 8T-, 15T-, 30T-, and 40T-period responses are observed, which shows the appearance of subharmonics. With the further decrease of f pre , chaotic behavior at ω = 0.23 is clearly visible in Figure 14d.

Conclusions
The surface topography has a significant effect on the vibration characteristics and stability of the contact interface system. In this article, an improved model of contact stiffness is developed to supply a study foundation. Furthermore, based on the contact model, a single-freedom dynamic model of the joint parts with the consideration of fractal surface has been established. In addition, the influences of excitation force and key design parameters of the machined bolt joint interface are discussed by amplitude-frequency curves, bifurcation diagrams, 3-D frequency spectrums, time histories, phase plane diagrams, and Poincare sections. Theconclusions are as follows: 1. The stiffness of the system increases nonlinearly with the deflection, which shows the hardening type nonlinearity in the positive direction and the softening type nonlinearity in the opposite direction. However, the system performs softening behavior under harmonic excitation. 2. The surface topography is controlled by fractal parameters D and G. With the increase of D or the decrease of G, the stiffness increases if the displacement is positive, and the bolt joint becomes relatively unstable. 3. The system is sensitive and unstable when it has a higher excitation force. When the bolt pre-tightening is controlled and the interface has not been separated, the joint parts present a poor stability with a smaller bolt pre-tightening force. Moreover, as the apparent contact area increases or the number of bolts decreases, subharmonics appear.

Conclusions
The surface topography has a significant effect on the vibration characteristics and stability of the contact interface system. In this article, an improved model of contact stiffness is developed to supply a study foundation. Furthermore, based on the contact model, a single-freedom dynamic model of the joint parts with the consideration of fractal surface has been established. In addition, the influences of excitation force and key design parameters of the machined bolt joint interface are discussed by amplitude-frequency curves, bifurcation diagrams, 3-D frequency spectrums, time histories, phase plane diagrams, and Poincare sections. Theconclusions are as follows: 1.
The stiffness of the system increases nonlinearly with the deflection, which shows the hardening type nonlinearity in the positive direction and the softening type nonlinearity in the opposite direction. However, the system performs softening behavior under harmonic excitation.

2.
The surface topography is controlled by fractal parameters D and G. With the increase of D or the decrease of G, the stiffness increases if the displacement is positive, and the bolt joint becomes relatively unstable. 3.
The system is sensitive and unstable when it has a higher excitation force. When the bolt pre-tightening is controlled and the interface has not been separated, the joint parts present a poor stability with a smaller bolt pre-tightening force. Moreover, as the apparent contact area increases or the number of bolts decreases, subharmonics appear.