Strojniški vestnik - Journal of Mechanical Engineering 62(2016)11, 665-677 © 2016 Journal of Mechanical Engineering. All rights reserved.
D0l:10.5545/sv-jme.2016.3410	Original Scientific Paper
Received for review: 2016-01-22 Received revised form: 2016-06-11 Accepted for publication: 2016-09-02
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
Yongsheng Zhao1* - Cheng Yang1 - Ligang Cai1 - Weimin Shi2 - Yi Hong1 1 Beijing University of Technology, Key Laboratory of Advanced Manufacturing Technology, China 2 Beijing University of Technology, College of Computer Science, China
Bolted connections are widely employed to fix the structural components, in which the bolted joint is one of the weakest parts and can significantly affect the dynamic characteristics of the machine tool. In this research, a stiffness and damping model based on the uneven surface contact pressure is presented for the bolted joint to accurately predict the dynamic characteristic of a bolted assembly. The normal and tangential stiffness and damping of the contact surface can be deduced based on the fractal contact theory. However, the contact surface pressure of bolted joint is unevenly distributed due to the influence of the concentrated force of multi-bolts. Therefore, the pressure of the contact surface is introduced to define the stiffness and damping of bolted joint. The assumption is that the contact surface is flat in the macro-scale. Then, we can obtain the pressure distribution of contact surface through the finite element (FE) method. The nonlinear relationship of stiffness, the damping of the bolted joint, and the pressure of contact surface can be obtained and assigned to the FE model based on the pressure distribution of the contact surface. An experimental set-up with a box-shaped specimen is designed for validating the proposed model. The equal pre-tightening force and bending moment effect case studies are provided to demonstrate the effectiveness of the model. The results show that the proposed model can be used to accurately predict the dynamic characteristic of the machine tool. Keywords: bolted joint, contact stiffness and damping, fractal contact theory, uneven pressure distribution, machine tool
Highlights
•	A stiffness and damping model based on the uneven surface contact pressure is presented for the bolted joint to accurately predict the dynamic characteristic of a bolted assembly.
•	The surface contact pressure is introduced to define the stiffness and damping of bolted joint, in which the uneven pressure distribution of contact surface can be obtained through the finite element (FE) method.
•	The nonlinear relationship of stiffness, the damping of the bolted joint, and the pressure of contact surface can be obtained and assigned to the FE model based on the pressure distribution of the contact surface.
•	An experimental set-up with a box-shaped specimen is designed for validating the proposed model.
•	The equal pre-tightening force and bending moment effect case studies are provided to demonstrate the effectiveness of the model.
0 INTRODUCTION
It is well known that the dynamic behaviour of bolted joints plays a significant role in determining the machining stability and precision of a machine tool. The dynamic behaviour of bolted joint can be affected by the geometry and machining precision of the surface, the property of material, external load, and pre-tightening force of bolts [1]. All of those factors and the complicated contact mechanism make it difficult to build an accurate model of bolted joints. The bolted joint also shows a strong nonlinear property, and it depicts with difficulty the relationship of the dynamic characteristic of bolted joints and influencing factors. The displacement measurement techniques can be used to obtain the frictional properties of bolted joint, such as LVDT Method [2], extensometer method [3], digital image correlation method [4] and [5], laser interferometer [6] or ultrasound method [7]. However, the contact surface of a bolted joint is difficult to measure directly by using
those experimental methods. Recently, the receptance coupling substructure analysis (RCSA) method was adopted to identify the stiffness and damping of bolted joints [8] to [10]. However, the RCSA method is easily affected by the measuring noise and environment. The identified results could not be applied to the bolted assembly with different technological parameters.
The micro-contact theory is another method to obtain the contact stiffness and damping of bolted joints, which make it possible to accurately model bolted joints and to optimize bolted assemblies. The micro-contact mechanics of bolted joints is based on the geometry topography of surface and Hertz contact theory. There are three methods in geometry topography of surface: the statistics method [11], the fractal geometry method [12], and the reconstruction of experimental data method [13]. The statistics method can be affected by the resolution and sampling length of the measuring instrument, which do not adequately reflect all features of the rough surface. However, the actual surface is always multi-scaled
*Corr. Author's Address: Beijing University of Technology, Key Laboratory of Advanced Manufacturing Technology, China, yszhao@bjut.edu.cn
665
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
with non-uniform asperities. The reconstruction of the experimental data method is also limited by the testing instrument. In contrast, the fractal geometry method has the merit of independent resolution and sampling length of measuring instrument and is suitable for the machined surface [14]. Therefore, researchers have focused on studying the contact mechanics of surface based on the fractal contact theory.
For the normal and tangential contact stiffness of bolted joints, most researchers have studied the contact model of asperities based on the Hertz contact theory and obtained the contact stiffness by integrating each micro-contact in the contact surface. Majumdar and Bhushan [14] and [15] introduced scale-independent fractal parameters to characterize the contact surface topography, known as the Majumdar-Bhushan (MB) model. Zhang et al. [16] built up the fractal model for normal contact stiffness of machine joint surfaces based on the Hertz contact theory between a sphere and a plane. Jiang et al. [17] used three different machining methods (milling, grinding, and scraping) iron specimens to obtain the contact stiffness under different contact loads. Shi et al. [18] established the contact stiffness and micro parameters on the joint surface based on the Greenwood and Williamson theory. Sherif [19] found that the ratio between the normal and tangential contact stiffness is constant and is merely a function of Poisson's ratio of the contact surfaces. Fuadi et al. [20] estimated the tangential contact stiffness of the contact interface with controlled contact asperities based on an experimental method.
For the normal and tangential contact damping of bolted joints, the energy dissipation model of bolted joints is adopted to obtain the contact damping. Zhang and Ding [21] proposed a fractal model of normal damping for joint interfaces based on the modified MB fractal model and mechanism of normal contact damping dissipating energy. Bograd et al. [22] studied tangential damping parameters based on a generic isolated joint test bench. The influence of joint parameters, such as normal contact force and frequency dependence is examined. Li et al. [23] proposed the fractal prediction model of the tangential contact damping of joint surface based on the "solidgap-solid" contact model and the Hertz theory. Zhang et al. [24] also proposed tangential contact damping of joint interfaces by considering the mechanism of energy dissipation and the elastic-plastic deformation regime of joint plane interfaces. The studies show that the stiffness and damping model of contact surfaces have a close relationship with pre-tightening force, geometry, material and other influencing factors,
and could be modeled by using the fractal method. Nevertheless, it does not correspond well with the case of a larger contact area and sparse distribution of multi-bolts.
The main objective of this study is to obtain the normal and tangential stiffness, and the damping of the bolted joint with uneven pressure distribution. The contact surface is assumed as flat in macro-scale, and the finite element method is used to obtain the uneven pressure distribution. The relationship of stiffness, damping, and pressure of contact surface can be deduced based on the fractal contact theory. The finite element model of assembly with bolted joint can be established, in which the MATRIX27 element is introduced to define the bolted joint. An experimental set-up with two box-shaped specimens is designed for validating the presented model. By comparing the experimental and theoretical natural frequency and mode shape, the maximum error of presented model is less than 5.1 %. The equal pre-tightening force and bending moment effect cases are also considered to demonstrate the effectiveness of the presented model.
1 NORMAL AND TANGENTIAL STIFFNESS, DAMPING MODEL BASED ON THE CONTACT PRESSURE
The normal, tangential stiffness and damping model is widely adopted to depict the bolted joint. It is usually assumed to be equal for each stiffness and damping element in the overall contact surface during the process of modelling the bolted assembly [8] to [10]. However, this hypothesis is not suitable for the bolted joint because of the influence of the concentrated force of multi-bolts. The uneven pressure distribution makes the stiffness and damping of bolted joint different for each other in the contact surface. In this section, the contact pressure is assumed to be evenly distributed locally but unevenly distributed across the whole contact surface. The contact pressure of surface is used as a variable to deduce the stiffness, damping of the contact surface.
1.1 Normal Stiffness and Damping Model
Fractal geometry theory can describe profile features of a rough surface and thus provide a means of characterizing asperities of a large range of sizes. Previous studies [12] to [14], have shown that machined surface can be simulated by the Weierstrass-Mandelbrot (W-M) function as:
z( x) = G(D-1) £
f f\ n \
cos2ny x
„(2-D ) n
1 < D < 2, y > 1, (1)
0
666
Zhao, Y. - Yang, C. - Cai, L. - Shi, W. - Hong, Y.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
where, n is the frequency index, nmax is the upper limit of frequency, Z is the height of the machined surface, z is the lateral distance, y denotes the scaling parameter for determining the spectral density and self-affine property, G is the fractal roughness parameter where a higher value of G indicates a rougher surface. D denotes the fractal dimension of a surface profile, which determines the contribution of high and low-frequency components in the surface profile. The fractal parameters D and G can be obtained by the power spectral density method based on the measured data of surface.
Nominal flat surfaces are rough in micro-scale. When brought into contact, an asperity is plastically or elastically deformed depending on its size. For simplicity, the two contact rough surfaces are supposedly equivalent to an assumed rigid flat surface and a rough surface, as shown in Fig. 1a. The real contact area is only a very small part of the nominal contact area (< 1 %). It is assumed that the numerous circular asperities are sufficiently apart from each other, and the asperities interactions are neglected [24]. The statistical distribution of the truncated microcontact area a' can be described as [6]:
n(d) = -Dy a'LD 2a'-(2+D)/2 for 0 < d< dL, (2)
where, a' = 2a represents the truncated area of a micro-contact, and a is the real contact area of a micro-contact, a'L is the truncated area of the largest elastic micro-contact, which is shown in Fig. 1b, y describes the domain extension factor for microcontact size distribution and decided by the equation:
|y2-D)/2 _ (1+¥d'2)-(2-D)/d]/[(2-D) /D] = 1. (3)
a)
b)
Fig. 1. Contact mechanism; a) contact of a rough surface and a flat surface, and b) contact region of a single asperity
The critical truncated area demarcating the elastic and plastic deformation regimes can be expressed as [24]:
a' = 2a =2
2E
H
. n2/(D-1)
G2
(4)
where, H = 2.8 Y is the hardness of the softer material [25], Y is the yield strength of softer material. Asperities with a > a' are in the full elastic deformation regime, and the asperities are in plastic deformation if they satisfy a < a'.
Based on Hertz contact theory, the deformation 8 and the equivalent radius R for one micro-contact are determined by the W-M fractal function as [8]:
5 = 23-DG (d-1) (ln yf n(D-2)l2 (df^'2, (5)
r,'D> 2
R =-
24-DnD 2GD-1 (ln yf
(6)
For the cases of elastic or fully plastic asperity deformations, the relationship of contact force f f and truncated area a' for one micro-contact can be given by [16].
f =- E * R 2S2 = e 3
29/2-dE*
n(D-3)l2 (ln yf G(D-1)a'(3-Dy2, (7)
fP = Hd,
(8)
where the subscripts e and p denote elastic and full plastic	deformations, _	respectively.
E * = ((l _ V[2 )E1 + (l _ v22 )/E2) 1 is equivalent elastic modulus, El, E2, v1 and v2 are the elastic modulus and Poisson ratios of two surfaces, respectively.
The total elastic force FE, plastic force FP, and real contact area Ar can be obtained by integrating the asperities over the contact surface.
F„ =
J fen(d)dd
2(9-2 D)2 D
3(3 - 2D)(3-D)/2
E' (ln y) G1
I/2 g(D-1)
D D ( 3
x y/ 2 a'L 2
3-D
(D * 1.5), (9) 2EVn^ (lny)/ G2dLMn(D = 1.5)
1 3
3
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
667
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
'	un, (2-d)2
Fp = f; fpn(a')da' = q ¿^alf^2, (10)
Ar = | c a'n(a')da' +lJ 1 dn(a')da' =
n 1-D
LJ	2 tD/li rl-D/2 , fl-Di2 \
= 4_2D¥ 2aL (l + a (11)
The total force F of contact surface is composed of the elastic and plastic force. Based on the definitions of pressure in the contact surface, the mean pressure of contact surface in the nominal area Ar can be written as:
A
2(9-2D)2DE' (lnyf G(D-V~-3(3 - 2D)An((-D)/2
(2-D)
HDy 2
j —D —D
xl a',2 - a'2 1+
(2 - D )
I ill j ' 2E>4 (lny)) G2a'L4 lnl ^
(2-D)
t'LDl2a'c 2 , (D * 1.5) (12)
n4 A
I
2HDw4 ,I , , .
+-— aL 4 a'4, (D = 1.5)
A
The normal contact stiffness kn for each elastic micro-contact can be written as:
f=dma_= -a) (13)
n dS dS/da' 3y[2n(2 - D)
Since all asperities are plastic deformations when a < a', the normal contact stiffness should take kn = 0. The whole normal stiffness KN of contact surface can be given as:
Kn =| a knn(d)da = 4E' (3 - D)D
3yjln (l - D)(2 - D)
¥l-°12a'LDI2 (( - a'c(1-D)l2
(14)
The asperities of contact surface are elastic or plastic deformations. The elastic energy is generated due to the effect of stiffness in the elastic deformation region, and the plastic energy is generated due to the effect of damping in the plastic deformation region. It is plastic deformation for micro-contact when a'< a' ,
the plastic energy of one micro-contact is given by [21]:
wnP = fS fpdS = Ha'S.
(15)
The plastic strain energy in plastic contact region can be obtained as:
WNP = Wnpn(a')da' =
2 2-D TJT^r^D-1 2-D	J D-2 D
= 2 D V 2 (ln Y) n 2 aL^'r. (16)
When the truncated area a' of a micro-contact satisfies a' < a' < aL , the elastic energy of one microcontact can be written as:
= f fdS = — E'R2 S Jo e 15
1 5 2
(17)
The elastic energy of contact surface can be obtained by integrating Eq. (17), results in:
7ne =| a wnen(a')da' =
"-2D * 2( D1) 1-D
22 E DG2(D-1 y 2 In y ( ,
5	(aL
). (18)
1 5(5 - 3D)n2
The normal contact damping dissipation factor can be obtained as [21]:
15(5 - 3D) Hn2 at
V 13-d	- —
NE 22 E Gd-1(2 - D)(ln Y)2(a'L 2 - a
3
TD '2-D
2 n 1 5-3D	5-3D
(19)
- ' 2 )
The normal damping of contact surface can be written as:
cn =n„MN
15(5 - 3D)Hn2 a'2
22 E*Gd-1(2 - D)(ln Y)2(a'L 2 - a' 2 )
M if (3 D f-yl-D'2a'rD'2 (a™ - a'(1,D)'2), (20)
3^(1 -D)(2-D f L { L	c >
where, M is the mass of the structure.
1.2 Tangential Stiffness and Damping Model
When the contact surface suffers tangential force under a certain normal force, the deformation of asperities on the surface will occur in the shear direction. According to [25], the tangential deformation of the single asperity can be given as:
668
Zhao, Y. - Yang, C. - Cai, L. - Shi, W. - Hong, Y.
m
D
+
NP
13
1 5-3D
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
8 =
3/jp
16G ' r
1 -11 -
MP
2/3
(21)
where p, t is the normal and tangential load of a single asperity respectively, G is the equivalent shear modulus of the contact surfaces, 1/G' = (2-vl)/G1 + (2-v2yG2, where Gb G2 is he modulus of the two materials, respectively. ^ is the static friction factor. r is the contact radius of one micro-contact.
For one micro-contact, the ratio of shear force and normal force satisfies t/p = T /P, T is the tangential force acting on the contact surface and described as T=TbAr , where Tb is the shear strength of the softer material [26]. The tangential stiffness of one microcontact can be expressed as kt=dt/ ddt, then the total tangential stiffness can be written as follows:
K-r =
J 1 ktn(a ')da ' = 8G'D¥(2-D)I2 f1 _ T
V2n(1 _ D) ^ HP
13
xa,
(DI2
(ai
t(1_D)l2
t(1_D)l2

(22)
For one single micro-contact, a small relative motion will be caused when a tangential force is less than the limiting friction force, as depicted in [24]. The contact area of micro-contact can be divided into a "microslip" section and a "stick" section. The hysteretic damping of contact surface mainly depends on the microslip of asperities for the contact surface, which is the main reason for the energy dissipation of contact surface.
According to [24], the energy dissipation of one micro-contact per cycle in the vibration can be expressed as:
_ 9n1/2«VP2
10GA'
xU-11 -
5T
ß.P ) 6ßP
1 + 11 -
ßP
2/3 \
>. (23)
The energy dissipation of contact damping for the whole surface can be expressed as:
Wd =| a wdn(a')da'
=_9 (Inf y2 P2 (4 - 2D )2_
20GV(2-dV2 D (3 - D )af2 (a(2 + af-^2 )2
/MA - a'(3-D)/2 ) I j-| j I3 5T
yP) 6yP
1 + 11
2
T
yP )
(24)
The total energy input W is a combination of dissipation energy Wd and elastic strain energy We of the two contact surfaces. Then, the elastic strain energy of the two surfaces can be obtained by the difference of total energy storage and dissipation energy per cycle. The energy storage of one microcontact w caused by the tangential force can be expressed as [27]:
w = 2n
3ß2 a2 P2
16 A-G'r
^	rp 5	rp 2
- + -(1 - —)- - (1 - —)-5 5 ßP	ßP
. (25)
The total energy input W per cycle for the contact surface due to tangential force can be obtained as:
W _ f "Lwn(a')da'
J a'c
__3n(2nf y2 P2 (4 - 2 D )2_
32GV(2-D)2 D (3 - D ) af2 (2-D))2 + af-D)2 )

(2-D)2 - a'(3-D)2
3+2 it -
5 51 ßP
5
T A 3
-it -
2
T V
VP,
. (26)
The elastic strain energy of two surfaces We can be written as:
W = W - Wd.	(27)
According to the definition of damping dissipation factor, it is the ratio of the dissipated energy to total elastic strain energy. The tangential contact damping dissipation factor can be expressed as:
n =-
(28)
The tangential contact damping mainly depends on hysteretic damping, which is the result of energy dissipation in the contact surfaces. The total tangential contact damping can be written as [24]:
x
d
d
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
669
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
Fig. 2. Flowchart diagram for the stiffness and damping model of the bolted joint
CT = ntKT =
1
(5n/24 )(H2/ H1)-1
8G' Dw
(22-D )/2
V2n(1 - D)
1 T I
1--I a
. VP J
V/3
,D/2( a ,(1-D)/2 L ( L
t(1-D )/2
).(29)
where, Hi, H2 are two nonlinear functions of T and P, which can be expressed as:
H = 1 -
1 -
uP
5T
H 2 = 3 + 1 2 5 5
6uP
T
T ^ 3
uP
ßP
ßP
In this section, the local uniform pressure case is introduced to study the normal and tangential stiffness, damping of contact surfaces based on the fractal contact theory. However, the contact pressure of bolted joint is uneven distribution due to the effect of the concentrated force of multi-bolts. To obtain the contact pressure distribution of surface, the contact surface is assumed to be flat in macro-scale. Then, the nodes of the contact surface are selected, and the contact stress can be extracted from the general Postproc section of ANSYS software. Based on the nodal stress of contact surface, the largest truncated area aL can be calculated by using Eq. (12). The normal and tangential stiffness, damping of the contact surface can be obtained by using Eqs. (14), (20), (22) and (29) based on the obtained area a, . The
values of stiffness and damping are assigned to the MATRIX27 element, which is used to connect the node-to-node of two contact surfaces. Fig. 2 shows the flowchart diagram for the stiffness and damping model of the bolted joint.
2 EXPERIMENTAL SET-UP AND VALIDATION FOR THE BOLTED JOINT
In order to validate the presented model of the bolted joint, a box-shaped specimen is designed as shown in Fig. 3. The material of the specimen is nodular cast iron (material type QT600-3). The contact surfaces are machined by milling with roughness Ra = 1.6 fim. The material property of the specimen is listed in Table 1. A surface profilometer is used to measure the contact surface for obtaining the data of surface profile. The fractal dimension D and fractal roughness parameter G can be computed based on the power spectrum density method [17], the reticular cell counting method [28], the variation method [29], yardstick method [30] and structure-function method [31]. In this paper, the power spectrum density method is used to obtain the fractal dimension D and fractal roughness parameter G. The power spectral density P(a>) of W-M function can be calculated by the following equation as:
G 2(D-1)
P ((D) = G-D-5
v } 2ln y
(30)
where y is the scaling parameter for determining the spectral density and self-affine property y = 1.5, and rn is the spatial frequency. Taking the logarithm of Eq.
670
Zhao, Y. -
 Yang, C. - Cai, L. - Shi, W. - Hong, Y.
5


Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
(30), the linearized equation can be obtained, and its slope k and intercept b can be given by:
k = 2 D - 5,
b = 2(D -l)lg G - lg(2lny),
(31)
The power spectrum density of a specimen surface and its linear fitting function can be obtained by adopting the least square method, as shown in Fig. 3. Based on the slope and intercept of the fitting function, the fractal dimension and fractal roughness parameter can be obtained as D = 1.42, G = 1.21E-8 m.
Fig. 3. Power spectrum density of a specimen surface
Fig. 4. Plane dimension of box-shaped specimen
The assembly is composed of two specimens and connected by eight M16 bolts. The material of bolts is 45# steel. The assembly is suspended to simulate the free degree of a freedom state, as shown in Fig. 5. Two group sensors are placed in two specimens, respectively. Each group consists of nine piezoelectric acceleration sensors. The sensor type is PCB Mode 330B30. An impact hammer is adopted for excitement. The signals of the impact hammer and the sensors are collected and analysed in an LMS modal analyser. In the experiment, the key problem is how to ensure an accurate pre-tightening force for each bolt. The errors
of the pre-tightening force will affect the contact status and dynamic characteristic of the bolted joint. It is difficult to obtain the accurate pre-tightening force only using the torque wrench, due to the effect of friction and the coupling effect of multi-bolts. Therefore, the pre-calibrated tension bolts are used to provide the accurate pre-tightening force.
Fig. 5. Experimental set-up of bolted joint
Table 1. Material property of two specimens
Parameter name	Parameter value
Density p (kg/m3)	7800
Elastic modulus E (Pa)	1.5E11
Poisson modulus	0.28
Hardness H (Pa)	1.96E9
In order to verify the presented stiffness and damping model based on the surface contact pressure, a finite element model for the assembly bolted joint is established in ANSYS as shown in Fig. 6. The contact pressure distribution of the bolted joint can be obtained based on the static analysis method. The pre-tightening force can be loaded in the ® to ® for each bolt as shown in Fig. 6. The nodes of the contact surface are selected, and the contact pressure could be extracted from analysis results. Then, we can use the values of nodes pressure to compute the contact stiffness and damping in the normal and tangential directions as depicted in section 2. Fig. 7 shows the mesh of contact surface with 1019 nodes, and two surfaces have identical meshes so that the nodes are in one-to-one correspondence on two surfaces, which makes it possible to add the self-defined element between two nodes. Here, the MATRIX27 element of ANSYS software is introduced to define the normal and tangential stiffness, damping of bolted joint respectively. The MATRIX27 stiffness element for a pair of nodes could be expressed as follows,
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
671
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
Kxx	0	0	0	0	0	-Kxx	0	0	0	0	0"
0	Kyy	0	0	0	0	0	-Kyy	0	0	0	0
0	0	Kzz	0	0	0	0	0	-Kzz	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0
- Kxx	0	0	0	0	0	Kxx	0	0	0	0	0
0	-Kyy	0	0	0	0	0	Kyy	0	0	0	0
0	0	-Kzz	0	0	0	0	0	Kzz	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0
0	0	0	0	0	0	0	0	0	0	0	0
where Kxx, Kyy are the tangential stiffness along the contact surface, Kzz is the normal stiffness perpendicular to the contact surface. The MATRIX27 damping element has a similar expression.
Fig. 6. FE model of assembly with bolted joint
Fig. 7. Mesh of contact surface
3 RESULTS AND DISCUSSION
The dynamic characteristic of a bolted joint can be affected by the contact status of the surface. Usually, the initial status of a bolted joint is always different with its operating status. The main reason is that
the bolted parts of the machine tool are affected by the external loads. As a result, the contact status of bolted joint can be changed by those external loads. Therefore, two typical cases are introduced to verify the presented method in this paper: the bolted joint with the equal pre-tightening force for each bolt and the bolted joint with bending moment effect. Figs. 8 and 9 show the contact pressure distribution of bolted joint. The pressure of contact surface is an uneven symmetrical distribution for the case of equal pre-tightening force, as shown in Fig. 8. For example, the maximum pressure appears near the bolt hole with 7.98 MPa and the pressure gradually declines to 0.62 MPa in the middle of two bolt holes when the pre-tightening force of each bolt is 9 kN. Fig. 9 shows the pressure distribution of contact surface with bending moment effect. The bending moment acting on the assembly alters the stress of the contact surface. An equivalent model can be presented to represent the bending moment effect by adjusting the pre-tightening force of each bolt. For example, the pre-tightening force of bolts increase from Fb to Fb' for ® to © and decrease from Fb to Fb" for © to ®, where Fb is the pre-tightening force without the effect of bending moment, Fb and Fb" are the pre-tightening force with the effect of bending moment, respectively. The
Fig. 8. Contact pressure distribution with the identical pre-tightening force for each bolt
Fig. 9. Contact pressure distribution with bending moment effect
672
Zhao, Y. - Yang, C. - Cai, L. - Shi, W. - Hong, Y.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
equivalent bending moment acting on the assembly can be written as M= La/ 2 (Fb'- Fb")Nb, where La is the distance between two-row bolts and Nb = 3 is the number of bolts in each side of the assembly. The pressure distribution is more non-uniform, in which the pressure of one side of the contact surface is improved, and that of the other side is decreased. Compared with the pressure distribution of the contact surface with only an equal pre-tightening force effect for each bolt, the bending moment can significantly affect the pressure distribution of the contact surface.
The finite element model of assembly with bolted joint can be established by assigning the values of stiffness and damping between two surfaces based on the contact pressure. In order to verify the presented stiffness and damping model based on the surface contact pressure for the bolted joint, the frequency response analysis is necessary for the presented model. The hammer mode experiment is adopted to obtain the experimental mode. The points of A and B are chosen as the exciting point, and points C, and D points are chosen as the pick-up point in the normal and tangential directions, respectively. The tangential and normal displacements of bolted assembly can be obtained. The frequency responses in a normal and tangential directions at pre-tightening force 15 kN for each bolt are shown in Figs. 10 and 11. Comparing the predicted frequency response with the measured one, it is almost coincident in the frequency range.
100 200 300 400 500 600 700 800 900 1000 Frequency (Hz)
Fig. 10. Frequency response in normal direction
	-Measured frequency response	
	-------Predicted frequency response	
		' 11 A
		jy\ ¡Jy V
		
		vy
100 200 300 400 500 600 700 S00 900 1000 Frequency (Hz)
Fig. 11. Frequency response in tangential direction
Table 2 shows the numerical and experimental 1st
to 6th order mode shapes. The experimental mode
shapes are measured by a modal hammering impact
testing method and the average results of three times
for the existing test are adopted to reduce the random
error and improve the signal-to-noise ratio. In order to
ensure the same pre-tightening force for each bolt, the
numerical constant torque wrench is used to adjust the
pre-tightening force of bolts. The results show that the
mode shapes from the introduced method have good
agreement with those of the experiment. In order to
verify the effect of uneven pressure distribution in the
contact surface, the uniform pressure distribution case
is introduced as depicted in [11] The surface pressure
for the uniform pressure distribution case can be
N I
obtained by the equation of pu = ^ Pn A , where, Nu
is the number of bolts, Pn is the pre-tightening force for each bolt, A is the whole area of bolted joint.
The 1st to 6th order natural frequencies of structure can be obtained as shown in Table 3. For obtaining the accurate experimental results, the specimens are decomposed and assembled three times again in the experiment. Every time the maximum error of pre-tightening force is less than 0.1 kN for each bolt. Furthermore, comparing the natural frequencies of three experiments, the maximum error is 0.41 %, which indicates good agreement for the bolted joint. In Table 3, the experimental natural frequency is the mean value of the natural frequency of three experiments. For the uniform pressure case with 9 kN pre-tightening force of each bolt, the errors of 1st to 6th order natural frequencies of assembly are 5.47 %, -8.62 %, -2.43 %, -8.1 %, -6.7 % and 4.7 %, respectively. Corresponding the errors of 1st to 6th order natural frequencies of assembly are 1.09 %, -3.79 %, 1.74 %, -3.88 %, -2.66 % and 2.48 % for the non-uniform pressure case. Comparing the natural frequency of the uniform pressure case, those of non-uniform pressure case are much closer to the experimental results. The results show that the uneven pressure distribution of contact surface can seriously affect the accuracy analysis of the bolted assembly. The contact stiffness and damping of the bolted joint can be changed by improving bolted pre-tightening force, which can enlarge the natural frequency of assembly. With the hammering test, the first order natural frequency is 526.03 Hz and 550.83 Hz for the pre-tightening force 9 kN and 15 kN of each bolt, respectively. The dynamic characteristic of the bolted assembly can be significantly improved with the increasing of the pre-tightening force.
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
673
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
When the bolted assembly is loaded by the bending moment, the contact pressure distribution is altered, as shown in Table 3. The pre-tightening force for each bolt is 15 kN, and the bending moment is M = 4500 Nm. An equivalent model is adopted to represent the bending moment effect by adjusting the pre-tightening force of each bolt. The pre-tightening forces of the bolts are 20 kN for ® to ®, 15 kN for @ and © and 10 kN for © to ®, as shown in Fig. 6. In comparison with the case without a bending moment, the 1st to 6th order natural frequency are depressed, which shows that the contact stiffness of the bolted joint is weakened by the bending moment. The main reason is that the uneven distribution of contact pressure is intensified due to the effect of the bending moment. The results show that it is unable to reflect
the effect of bending moment and has much error for the uniform pressure case. However, it still maintains good consistency with the experimental results for the non-uniform stress case. The presented stiffness and damping method based on the contact surface pressure is more accurate in the modeling of the bolted joint, whether or not considering the influence of external load. The presented method can be used to accurately predict the behavior of the bolted assembly. It is also found that the first order natural frequency can be increased with the improvement of the pre-tightening force, as shown in Fig. 12. However, the first order natural frequency increases slowly when the pre-tightening force is larger than 24 kN. The main reason is that the contact pressure near the bolts can increase with the increasing of pre-tightening force;
Table 2. Comparison of theoretical and experimental mode shapes (total displacement)
674
Zhao, Y. - Yang, C. - Cai, L. - Shi, W. - Hong, Y.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
Table 3. Comparison of natural frequency for experimental result, uniform and non-uniform pressure adopting the presented method
Pre-tightening force for each bolt		fi	f2	f3	f4	f5	f6
	Experimental results	526.03	666.09	820.16	897.98	1334.41	1373.73
	Uniform pressure case	556.5	613.21	840.12	830.64	1250.63	1441.49
9 kN	Error of uniform pressure case [%]	5.47	-8.62	2.43	-8.1	-6.7	4.7
	Non-uniform pressure case	531.74	640.83	834.67	863.05	1299.8	1408.7
	Error of non-uniform pressure case [%]	1.09	-3.79	1.74	-3.88	-2.66	2.48
	Experimental results	550.83	668.11	871.12	903.26	1390.67	1414.80
	Uniform pressure case	588.61	619.16	835.20	838.19	1263.35	1467.82
15 kN	Error of uniform pressure case [%]	6.42	-7.9	-4.3	-7.76	-10.07	3.6
	Non-uniform pressure case	555.26	641.8	854.63	872.8	1319.8	1407.5
	Error of non-uniform pressure case [%]	0.8	-3.94	-1.89	-3.37	-5.1	-0.52
	Experimental results	547.49	667.59	834.22	900.36	1355.13	1380.96
15 kN (with bending moment effect M = 4500 N-m)	Uniform pressure case	588.61	619.16	835.20	838.19	1263.35	1467.82
	Error of uniform pressure case [%]	6.98	-7.82	0.12	-7.42	-7.26	5.92
	Non-uniform pressure case	551.81	637.53	836.72	878.42	1308.1	1416.8
	Error of non-uniform pressure case [%]	0.78	-4.5	0.29	-2.44	-3.59	2.53
nevertheless, the contact stress far away from the bolts is almost invariant for the bolted joint.
580
6 9 12 15 18 21 24 27 30 Pre-tightening force (kN)
Fig. 12. Effect of pre-tightening force on the natural frequency of assembly
4 CONCLUSIONS
In this paper, a normal and tangential stiffness damping model based on the contact surface pressure is presented for the bolted joint. The relationship of stiffness, the damping bolted joint, and the contact surface pressure can be deduced. An experimental set-up with a box-shaped specimen is designed for validating the proposed model. The equal pre-tightening force and bending moment effect case studies are provided to demonstrate the effectiveness of the model. The maximum error of natural frequency for the uniform pressure case is twice that of nonuniform pressure cases. The stiffness of bolted joint is weakened as the pressure distribution is more uneven,
caused by the bending moment. The results show that the presented model is effective in the modeling of bolted joints whether having the influence of external load or not, which can meet the requirement of accurately predicting the behaviour of the bolted assembly for the machine tool.
5 ACKNOWLEDGEMENTS
This work was supported by Beijing Municipal Natural Science Foundation (CN) (No. 3132004), National Natural Science Foundation of China (No. 51375025) and National Science and Technology Major Project of China (No. 2013ZX04013-011).
6 NOMENCLATURES
a	real contact area of one micro-contact [m2]
A0	nominal area [m2]
Ar	total real contact area [m2]
C	total normal damping
D	fractal dimension
E	equivalent elastic modulus [Pa]
Ei,	E2 elastic modulus of two surfaces [Pa]
f	normal load of one micro-contact [N]
F	total normal load of contact surface [N]
G	fractal roughness parameter [m]
G'	equivalent shear modulus [Pa]
G1, G2 shear modulus of two surfaces [Pa]
H	hardness of the soft material [MPa]
k	stiffness of one micro-contact [N/m]
K	total stiffness for contact surface [N/m]
L	sampling length, [m]
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
675
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
M	mass of the structure [kg]
n	frequency index nmax upper limit of frequency index
p	normal load of one micro-contact [N]
P	normal force of contact surface[N]
r	radius of the real contact region [m]
R	curvature radius [m]
t	tangential load of one micro-contact [N]
T	shear force of contact surface [N]
w	energy storage of one micro-contact [J]
wd	dissipated energy of one micro-contact [J]
W	energy storage of contact surface [J]
Wd	dissipated energy of contact surface [J]
x	lateral distance [m]
Y	yield stress of the soft material [Pa] z	height of machined surface [m]
8	deformation of one micro-contact [m]
Y	scaling parameter of the spectral density n	contact damping dissipation factor
vj, v2 Poisson ratio
^	static friction coefficient
am	mean pressure of contact surface [Pa]
Tb	shear strength of the soft material [Pa]
p	expand coefficient
Subscripts or superscripts: ' truncated section of one micro-contact c critical parameter demarcating the elastic and
plastic regimes e parameter of one micro-contact in the elastic regime
E parameter of contact surface in the elastic regime L parameter for the largest asperity n normal parameter of one micro-contact N normal parameter of contact surface p parameter of one micro-contact in the plastic regime
P parameter of contact surface in the plastic regime t tangential parameter of one micro-contact T tangential parameter of contact surface
8 REFERENCES
[1]	Mao, K.M., Li, B., Wu, J., Shao, X.Y. (2010). Stiffness influential factor-based dynamic modeling and its parameter identification method of fixed joints in machine tool. International Journal of Machine Tools and Manufacture, vol. 50, no. 2, p. 156-164, DOI:10.1016/j.ijmachtools.2009.10.017.
[2]	Mohd Tobi, A.L., Ding, J., Bandak, G., Leen, S.B., Shipway, P.H. (2009). A study on the interaction between fretting wear and cyclic plasticity for Ti-6Al-4V. Wear, vol. 267, no. 1-4, p. 270282, DOI:10.1016/j.wear.2008.12.039.
[3]	Fridrici, V., Fouvry, S., Kapsa, Ph. (2003). Fretting wear behavior of a Cu-Ni-In plasma coating. Surface and Coatings
Technology, vol. 163-164, p. 429-434, DOI:10.1016/S0257-8972(02)00639-4.
[4]	Kartal, M.E., Mulvihill, D.M., Nowell, D., Hills, D.A. (2011). Determination of the frictional properties of titanium and nickel alloys using the digital image correlation method. Experimental Mechanics, vol. 51, no. 3, p. 359-371. DOI:10.1007/s11340-010-9366-y.
[5]	Kartal, M.E., Mulvihill, D.M., Nowell, D., Hills, D.A. (2011). Measurements of pressure and area dependent tangential contact stiffness between rough surfaces using digital image correlation. Tribology International, vol. 44, no. 10, p. 11881198, DOI:10.1016/j.triboint.2011.05.025.
[6]	Filippi, S., Akay, A., Gola, M.M. (2004). Measurement of tangential contact hysteresis during microslip. Journal of Tribology, vol. 126, no. 3, p. 482-489, DOI:10.1115/1.1692030.
[7]	Mulvihill, D.M., Brunskill, H., Kartal, M.E., Dwyer-Joyce, R.S., Nowell, D. (2013). A comparison of contact stiffness measurements obtained by the digital image correlation and ultrasound techniques. Experimental Mechanics, vol. 53, no. 7, p. 1245-1263, DOI:10.1007/s11340-013-9718-5.
[8]	Kashani, H., Nobari, A.S. (2010). Identification of dynamic characteristics of nonlinear joint based on the optimum equivalent linear frequency response function. Journal of Sound and Vibration, vol. 329, no. 9, p. 1460-1479, DOI:10.1016/j.jsv.2009.11.007.
[9]	Celic, D., Boltezar, M. (2009). The influence of the coordinate reduction on the identification of the joint dynamic properties. Mechanical Systems and Signal Processing, vol. 23, no. 4, p. 1260-1271, DOI:10.1016/j.ymssp.2008.11.002.
[10]	Wang, L.Y., Yin, G.G., Zhang, J.F. (2006). Joint identification of plant rational models and noise distribution function using binary-valued observation. Automatica, vol. 42, no. 4, p. 535547, DOI:10.1016/j.automatica.2005.12.004.
[11]	Shi, X., Polycarpou, A.A. (2005). Measurement and modeling of normal contact stiffness and contact damping at the meso scale. Journal of Vibration and Acoustics, vol. 127, no.1, p. 5260, DOI:10.1115/1.1857920.
[12]	Ciulli, E., Ferreira, L.A., Pugliese, G., Tavares, S.M.O. (2008). Rough contacts between actual engineering surfaces: Part I. Simple models for roughness description. Wear, vol. 264, no. 11-12, p. 1105-1115, DOI:10.1016/j.wear.2007.08.024.
[13]	Komvopoulos, K., Ye, N. (2000). Three-dimensional contact analysis of elastic-plastic layered media with fractal surface topographies. Journal of tribology, vol.123, no.1, p.632-640, DOI:10.1115/1.1327583.
[14]	Majumdar, A., Bhushan, B. (1990). Role of fractal geometry in roughness characterization and contact mechanics of surfaces. Journal of Tribology, vol. 112, no. 2, p. 205-216, DOI:10.1115/1.2920243.
[15]	Majumdar, A., Bhushan, B. (1991). Fractal model of elastic-plastic contact between rough surfaces. Journal of Tribology, vol. 113, no. 1, p. 1-11, DOI:10.1115/1.2920588.
[16]	Zhang, X.L., Huang, Y.M., Han, Y. (2000). Fractal model of the normal contact stiffness of machine joint surfaces based on the fractal contact theory. China Mechanical Engineering, vol. 11, no. 7, p. 727-729. (in Chinese)
676
Zhao, Y. - Yang, C. - Cai, L. - Shi, W. - Hong, Y.
Strojniski vestnik - Journal of Mechanical Engineering 62(2016)11, 556-677
[17]	Jiang, S.Y., Zheng, Y, Zhu, H. (2009). A contact stiffness model of machined plane joint based on fractal theory. Journal of Tribology, vol. 132, no. 1, p. 1-6, D0l:10.1115/1.4000305.
[18]	Shi, J.P., Ma, K., Liu, Z.Q. (2012). Normal contact stiffness on unit area of a mechanical joint surface considering perfectly elastic elliptical asperities. Journal of Tribology, vol. 134, no. 3, p. 1-6, D0I:10.1115/1.4006924.
[19]	Sherif, H.A., Kossa, S.S. (1991). Relationship between normal and tangential contact stiffness of nominally flat surfaces. Wear, vol. 151, no. 1, p. 49-62, D0I:10.1016/0043-1648(91)90345-U.
[20]	Fuadi, Z., Takagi, T., Miki, H., Adachi, K. (2013). An experimental method for tangential contact stiffness evaluation of contact interfaces with controlled contact asperities. Proceedings of the Institution of Mechanical Engineers, Part J: Journal of Engineering Tribology, vol. 227, no. 10, p. 1117-1128, D0I:10.1177/1350650113481293.
[21]	Zhang, X.L., Ding, H.Q. (2013). Normal contact damping and dissipation factor model of joint interfaces based on fractal theory. Chinese Transactions of the Chinese Society for Agricultural Machinery, vol. 44, no. 6, p. 287-294, D0I:10.6041/j.issn.1000-1298.2013.06.050.
[22]	Bograd, S., Schmidt, A., Gaul, L. (2008). Joint damping prediction by thin layer elements. Proceedings of the IMAC 26th. Society of Experimental Mechanics Inc. Bethel.
[23]	Li, X.P., Wang, W., Zhao M.Q., Nie, H.F., Wen, B.C. (2012). Fractal prediction model for tangential contact damping of joint surface considering friction factors and its simulation. Journal of Mechanical Engineering, vol. 23, no. 12, p. 46-50, D0I:10.3901/jme.2012.23.046. (in Chinese)
[24]	Zhang, X.L., Wang, N.S., Lan, G.S., Wen, S.H., Chen, Y.H. (2013). Tangential damping and its dissipation factor models of joint interfaces based on fractal theory with simulations. Journal of Tribology, vol. 136, no. 1, p. 1-10, D0I:10.1115/1.4025548.
[25]	Kogut, L., Etsion, I. (2002). Elastic-plastic contact analysis of a sphere and a rigid flat. Journal of Applied Mechanics, vol. 69, no. 5, p. 657-662, D0I:10.1115/1.1490373.
[26]	Johnson, K.L. (1985). Contact Mechanics. Cambridge University Press, Cambridge, D0I:10.1017/ CB09781139171731.
[27]	Tian, H.L., Zhao, C.H., Fang, Z.F. Zhu, D.L., Li, X., Mao, K.M. (2013). Improved model of tangential stiffness for joint interface using anisotropic fractal theory. Transactions of the Chinese Society for Agricultural Machinery, vol. 44, no. 3, p. 257-266, D0I:10.6041/j.issn.1000-1298.2013.03.046.
[28]	Goerke, D., Willner, K. (2008). Normal contact of fractal surfaces-Experimental and numerical investigations. Wear, vol. 264, no. 7, p. 589-598, D0I:10.1016/j.wear.2007.05.004.
[29]	Boneh, Y., Sagy, A., Reches, Z. (2013). Frictional strength and wear-rate of carbonate faults during high-velocity, steady-state sliding. Earth and Planetary Science Letters, vol. 381, p. 127137, D0I:10.1016/j.epsl.2013.08.050.
[30]	Ai, T., Zhang, R., Zhou, H.W., Pei, J.L. (2014). Box-counting methods to directly estimate the fractal dimension of a rock surface. Applied Surface Science, vol. 314, p. 610-621, D0I:10.1016/j.apsusc.2014.06.152.
[31]	Wang, T., Wang, L., Zheng, D., Zhao, X., Gu, L. (2015). Numerical simulation method of rough surfaces based on random switching system. Journal of Tribology, vol. 137, no. 2, p. 021403, D0I:10.1115/1.4029644.
Stiffness and Damping Model of Bolted Joints with Uneven Surface Contact Pressure Distribution
677