Dynamic modeling method of the bolted joint with uneven distribution of joint surface pressure

Li, Shichao; Gao, Hongli; Liu, Qi; Liu, Bokai

doi:10.1063/1.5021823

The dynamic characteristics of the bolted joints have a significant influence on the dynamic characteristics of the machine tool. Therefore, establishing a reasonable bolted joint dynamics model is helpful to improve the accuracy of machine tool dynamics model. Because the pressure distribution on the joint surface is uneven under the concentrated force of bolts, a dynamic modeling method based on the uneven pressure distribution of the joint surface is presented in this paper to improve the dynamic modeling accuracy of the machine tool. The analytic formulas between the normal, tangential stiffness per unit area and the surface pressure on the joint surface can be deduced based on the Hertz contact theory, and the pressure distribution on the joint surface can be obtained by the finite element software. Futhermore, the normal and tangential stiffness distribution on the joint surface can be obtained by the analytic formula and the pressure distribution on the joint surface, and assigning it into the finite element model of the joint. Qualitatively compared the theoretical mode shapes and the experimental mode shapes, as well as quantitatively compared the theoretical modal frequencies and the experimental modal frequencies. The comparison results show that the relative error between the first four-order theoretical modal frequencies and the first four-order experimental modal frequencies is 0.2% to 4.2%. Besides, the first four-order theoretical mode shapes and the first four-order experimental mode shapes are similar and one-to-one correspondence. Therefore, the validity of the theoretical model is verified. The dynamic modeling method proposed in this paper can provide a theoretical basis for the accurate dynamic modeling of the bolted joint in machine tools.

The dynamic characteristics of machine tool influence the machining accuracy, the machine performance and the machine life significantly. Therefore, establishing an accurate dynamic model of machine tool is very important to evaluate and optimize the dynamic characteristics of the machine tool. Unlike a single component, a machine tool is an integrated system that composed of many components with a lot of joints. The joint is essentially a mechanical joint surface, which is the surface that the parts contact with each other. Various types of joints in the machine tool have important influence on the dynamic characteristics of the whole machine. Research shows that the stiffness of the joints is 60% ∼ 80% of the whole machine stiffness, and more than 90% of the damping of the whole machine originate from the joints.^1–3 Therefore, it is necessary to consider the influence of the dynamic characteristics of the joints on the dynamic characteristics of the whole machine during the dynamic modeling of the machine tool, and establishing a reasonable dynamic model of the joints will help to improve the accuracy of the whole machine dynamics model. However, there are many factors influence the dynamic characteristics of the joint (such as material properties, joint surfaces topography, preload, and whether there is medium between the joint surfaces) and the contact mechanism is complex. Therefore, the dynamic modeling of the joints has always been one of the difficulties in machine tool dynamics modeling. Through the investigation found that in various types of joints, the bolted joints are widely found in various types of machine tools. Therefore, this paper researches on the dynamic modeling method of the bolted joints.

In recent years, many scholars at home and abroad have done a lot of research on how to establish an accurate dynamics model of the joint, and obtained a lot of valuable research results. Greenwood and Williamson proposed the probability-statistical contact model of two contact surfaces, that is the GW model.⁴ The model considers the height distribution of the surface topography as a random variable. And through combining the Hertz contact theory and the statistical knowledge, the relationship between the pressure and deformation of the joint surface can be obtained. The GW model has five assumptions: (1) the rough surface is isotropy, (2) the asperity is sphere near the peak, (3) all asperities have the same radius but their heights distributed randomly, (4) the asperities are far away from each other, and there is no interaction between them, (5) the rough surfaces are not deformed, and only the asperity deform during contact. Because of the consideration of the surface topography parameters, a satisfactory explanation for the classical friction law can be gotten when the surface height obeys the one-dimensional normal distribution. So that the GW model is still accepted by the majority of researchers. Because the surfaces of many components have fractal properties, Majumdar and Bhushan⁵ proposed a rough surface elastic-plastic contact fractal model with scale independence which is based on the Weierstrass-Mandelbrot fractal function. The model describes the relationship between the normal force and normal deformation of the joint. In 1979, Yoshimura⁶ made a research on the joints of machine tools. He found that the dynamic force of a joint was a combination of six forces in different forms, that is the six DOFs generalized forces. An equivalent dynamic model of six DOFs between joint surfaces was established by simplifying the joint to a node, of which six DOFs are independent with each other. Furthermore, Youshimura derived the mathematical formulas for calculating the equivalent stiffness and equivalent damping in each DOF model. S. Taylor et al.^7–10 described the dynamic characteristics of the joint by establishing a series of spring-damping elements between the joint surfaces. The stiffness and damping of the spring-damping elements can be obtained by means of analytical calculations or parameter identification. Because the method is simple and practical, it is widely adopted by many scholars and engineers. Zhang Jie et al.¹¹ thought that the dynamic properties of the joints under various conditions (including different materials, surface roughness, pre-tightening force, medium, etc.) can be equivalent substituted with the dynamic characteristics of the ideal joint surface which have the same mechanical properties in any position. Based on this idea, a finite element is proposed to establish the ideal joint surface dynamic model. Meanwhile, according to the theory of finite element, the stiffness matrix and damping matrix of the element are deduced, and the parameters of which are obtained by parameter identification. H Tian et al.¹² simulated the dynamic characteristics of the joint surface as a layer of virtual material that is fixedly connected to the components of the joint, and the elastic modulus, density and poisson’s ratio of the virtual material are deduced based on the fractal theory and the Hertz contact theory. For the bolted joints, K Mao et al.¹³ developed a new element and it is used for joint modeling. The element has 8 nodes and 24 degrees of freedom. Then the stiffness matrix and damping matrix of the element are deduced, and the parameters of the matrices are obtained by parameter identification.

The bolted joints has a phenomenon that the pressure distribution of the joint surface is uneven because of the preload force of the bolts. Besides, research shows that the joint surface pressure is one of the important factors that determine the contact stiffness and damping of the bolted joints.^14–16 However, it can be found through the analysis of the above papers that these methods do not consider the uneven pressure distribution of the joint surface during the modeling process. So considering the joint surface uneven pressure during the dynamics modeling process of the bolted joints can effectively improve the modeling precision. Based on the phenomenon that the joint surface pressure distribution is uneven, this paper presents a general method to build bolted joints finite element model according to the joint surface uneven pressure distribution, and the validity of the modeling method is verified by comparison with the modal test data and the result of other modeling methods.

I. A new dynamic model of the bolted joints in machine tools

A. Description of the research object

After investigating various types of CNC lathes, an assembly specimen composed of two blocks is designed and manufactured. The blocks are connected by four M8 bolts, as depicted in FIG. 1a. The unit for the dimensions given in FIG. 1b is mm. The upper block and the lower block are connected by bolts, so the contact between the upper block and the lower block is a typical bolted joint. The material of upper block and the lower block is 45# steel which correspond to the material property are given in TABLE I. In this paper, the bolted joint in the assembly specimen is used as the research object to study the dynamic modeling method of the bolted joints.

FIG. 1.

View large Download slide

Plane dimension of the assembly specimen.

TABLE I.

Material property of assembly specimen.

Material name	Elastic modulus E(Gpa)	Poisson modulus σ	Density ρ(Kg/m³)
45# steel	210	0.31	7850

B. Dynamic modeling of the bolted joints

Taking the assembly specimen as an example, which is shown in FIG. 1. The dynamic modeling process of the bolted joints is presented as follows:

Step 1.
Meshing the assembly specimen in the ANSYS.
Step 2.
Based on the principle of the shortest distance, finding the closest grid node on the joint surface of the upper block to the reference point which is one of the grid node on the joint surface of the lower block. Then establishing a Matrix27 user-defined stiffness element of the ANSYS software between the two nodes, and repeat this process until all grid nodes on the joint surface of the lower block are established the Matrix27 stiffness element with the nodes on the joint surface of the upper block. In this way, the dynamic behaviour of the bolted joint is described.
Step 3.
Based on the Hertz contact theory, the functional relationships between the joint surface pressure and the normal or tangential stiffness per unit area on the joint surface are derived. Futhermore, the unknown parameters in the mathematical expression are identified based on modal test data and parameter identification method.
Step 4.
Applying the bolt preload to the finite element model of the assembly specimen, and making static analysis. After static analysis, the pressure distribution on the joint surface can be extracted from the general Postproc section of ANSYS software. Based on the functional relationships between the joint surface pressure and the normal or tangential stiffness per unit area, the normal and tangential stiffness distribution of the joint surface can be calculated. Furthermore, the values of normal and tangential stiffness are assigned to the Matrix27 stiffness elements on the joint surface, and the joint dynamics model is established at this point.
Step 5.
Compositing the finite element model of the bolted joint, the upper block and the lower block into the finite element model of whole machine.

The process of dynamic modeling of bolted joints is shown in FIG. 2.

FIG. 2.

View large Download slide

Flow chart of dynamic modeling of bolted joints.

1. Mesh generation

Finite element meshing of the assembly specimen was implemented by ANSYS software, which has strong meshing and numerical calculation capabilities. There are some curved surface such as holes in the assembly specimen. Therefore, in order to simulate the geometric characteristics and the stress concentration around the holes better, the Solid187 element in ANSYS software is chosen to mesh the assembly specimens. The Solid187 element is a tetrahedral element with 10 nodes, and it has a good geometric fitting property for the structure with curved surface. Futhermore, each node of the Solid187 element has three degrees of freedom. And it has two order displacement characteristics, which can simulate the stress variation inside the element and meets the accuracy requirements of engineering analysis. The assembly specimen after meshing is shown in FIG. 3.

FIG. 3.

View large Download slide

FE model of assembly specimen.

2. Establishing the Matrix27 stiffness elements on the bolted joints

Nominal flat surfaces are rough in micro-scale, as depicted in FIG. 4. The joint surface is not an absolutely smooth plane, but an uneven surface. When the components are connected together by threaded fasteners, many micro-contact pairs are formed by the bulges on the joint surfaces. The micro-contact pairs play important roles of transfering force and movement between two components. Moreover, they are random, even and discrete distribution on the joint surface. In the finite element model, the transfer of force and displacement between joints is achieved by grid nodes. The grid nodes are also random, even, discrete distribution on the joint surface, And this characteristic satisfies the distribution characteristics of micro-contact pairs on the joint surface. Therefore, it is reasonable to use finite element model to describe the force transfer characteristics of the joint, and the force transfer characteristics of the joint finite element model will approach the true force transfer characteristics of the joint with the increase of mesh density of the joint surface. Because the joint shows the ability to store and consume energy, each micro-contact pair can be described by a two-nodes element with stiffness characteristics to simulate the stiffness characteristics of the joint. For the bolted joints, the joint surfaces are tightly connected under the action of fasteners. Besides, the nonlinear behavior is restrained and not obvious, and the bolted joints shows linear characteristics. So the linear stiffness element can be established between the grid nodes of the bolted joints to describe the force transfer characteristics of the bolted joints. However, the grid nodes on the lower block joint surface are not necessarily coincide with the grid nodes on the upper block joint surface. The Matrix27 stiffness element is a linear stiffness element between two nodes with 6 degrees of freedom (including 3 horizontal degrees of freedom and 3 rotational degrees of freedom). It can describe the force transfer characteristics between two grid nodes with any distance, and it is chosen to describe the force transfer characteristics of the bolted joints in this paper.

FIG. 4.

View large Download slide

Micro-schematic of joint.

Taking the assembly specimen as an example, which is shown in FIG. 1. All grid node numbers and coordinates of the lower block joint surface are extracted from the general Postproc section of ANSYS software and represented as set A₁, as well as the data of the upper block is represented as set A₂. Selecting any grid node number in the set A₁ and figure out the distances between the grid node number and all grid node numbers in the set A₂ respectively. The calculation formula is shown in Eq. (1). Then establishing a Matrix27 user-defined stiffness element between the nearest two grid nodes, and repeat this process until all grid nodes in the set A₁ are established the Matrix27 stiffness element with the grid nodes in the set A₂.

L = \sqrt{{(x_{n}^{'} - x_{N})}^{2} + {(y_{n}^{'} - y_{N})}^{2} + {(z_{n}^{'} - z_{N})}^{2}}

(1)

where L is the distance between the two grid nodes. $x_{n}^{'}$ ⁠, $y_{n}^{'}$ and $z_{n}^{'}$ are the coordinate values of the grid node numbered n on the joint surface of the upper block. x_N, y_N and z_N are the coordinate values of the grid node numbered N on the joint surface of the lower block. After the above work, the establishment of the Matrix27 stiffness elements between the joint surfaces are finished. The positions of the nodes in all of the Matrix27 stiffness elements in the bolted joint is shown in FIG. 5, in which the circles represent the positions of the nodes that belong to the set A₁ in the Matrix27 stiffness elements and the red points represent the positions of the nodes that belong to the set A₂ in the Matrix27 stiffness elements. Statistically, there are 676 Matrix27 stiffness elements in the bolted joint.

FIG. 5.

View large Download slide

Nodes distribution map of the joint surface.

In this study, the work of establishing Matrix27 elements in the joint is completed by MATLAB and ANSYS. The detailed process is as follows:

Writing algorithms by MATLAB to calculate the distances from each node on one of the joint surfaces to all the nodes on the other joint surface respectively, and the calculation formula is shown in Eq. (1). Then marking the smallest distances as L, and judge whether L is smaller than ε. If the condition is satisfied, the two nodes are determined as a node pair. Therefore, establishing the Matrix27 stiffness element at each node pair in the joint based on the secondary-development programs by APDL language in ANSYS. Through the above work, all the Matrix27 stiffness elements are established.

3. Calculation of stiffness and damping coefficient

The Matrix27 user-defined stiffness element can be mathematically expressed as

[\begin{matrix} k_{xx} & 0 & 0 & 0 & 0 & 0 & - k_{xx} & 0 & 0 & 0 & 0 & 0 \\ 0 & k_{y y} & 0 & 0 & 0 & 0 & 0 & - k_{y y} & 0 & 0 & 0 & 0 \\ 0 & 0 & k_{z z} & 0 & 0 & 0 & 0 & 0 & - k_{z z} & 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 \\ - k_{xx} & 0 & 0 & 0 & 0 & 0 & 0 & k_{xx} & 0 & 0 & 0 & 0 \\ 0 & - k_{y y} & 0 & 0 & 0 & 0 & 0 & 0 & k_{y y} & 0 & 0 & 0 \\ 0 & 0 & - k_{z z} & 0 & 0 & 0 & 0 & 0 & 0 & k_{z z} & 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 \end{matrix}]

(2)

where k_xx, k_yy and k_zz are the stiffness between the two nodes in the x, y, and z directions respectively. And the stiffness are obtained by the functional relationships between the joint surface pressure and the normal or tangential stiffness per unit area on the joint surface. The rotational degrees of freedom in the Matrix27 element are neglected in this paper.

Research suggests that the stiffness of the bolted joint can be affected by the property of material, the size of joint surface pressure and the medium between two mating surfaces (such as lubricating oil, lubricating grease, etc.), especially the joint surface pressure has a significant effect on the stiffness of the bolted joints. In the bolted joints, the surface pressure of the joint surface is unevenly distribution due to the effect of the concentrated force of multi-bolts, so the stiffness on the joint surface is also unevenly distribution. According to the Hertz contact theory, the functional relationships of the normal and tangential stiffness per unit area and the surface pressure in the joint can be deduced.

a. Relationship between the joint surface pressure and the normal stiffness per unit area.

In the actual joint, the micro-contact pairs are formed by the contact of the bulges on the different joint surfaces. The microstructure of the bulges in the actual joint surface is usually an oval. The bulge can be considered as a sphere because the size of contact area on the oval is much smaller than its curvature radius.¹² The contact of two joint surfaces can be regarded as that of many spheres with different heights, and the contact of two spherical micro-contact pairs is shown in FIG. 6.¹²

FIG. 6.

View large Download slide

Elastic contact of micro-contact pairs.

In FIG. 6b, two spherical micro-contact pairs are flatten under the force F acting on the point O1 in the direction O1-O2. The distance between the spherical center O1 and O2 decreases the local interference w and it results in a contact circle with the radius r. The points A and B superpose a point C. The parameters R1 and R2 are the radii of micro-contact pairs 1 and 2 respectively.

According to the Hertz contact theory,¹⁷ the relationship between the contact force and the resulting deformation of a micro-contact pair can be described as follows:

F = \frac{4}{3} Ē \sqrt{R} w^{1.5}

(3)

where $Ē$ is the equivalent elastic modulus of the micro-contact pair, R is the equivalent radius of the micro-contact pair, w is the deformation of the micro-contact pair, F is the external force applied to the micro-contact pair. Futhermore, the $Ē$ and the R can be expressed as follows:

\frac{1}{Ē} = \frac{1 - μ_{1}^{2}}{E_{1}} + \frac{1 - μ_{2}^{2}}{E_{2}}

(4)

R = \frac{R_{1} R_{2}}{R_{1} + R_{2}}

(5)

Where E₁, μ₁ and R₁ represent the elastic modulus, the Poisson’s ratio and the radius of curvature of the micro-contact 1, and E₂, μ₂ and R₂ represent the elastic modulus, the Poisson’s ratio and the radius of curvature of the micro-contact 2.

According to Eq. (3), the normal stiffness formula of the micro-contact pair (normal force-deformation formula) can be expressed as follows:

K_{n} = \frac{d F}{d w} = 2 Ē \sqrt{R} w^{\frac{1}{2}}

(6)

where K_n is the normal stiffness of the micro-contact pair. According to Eq. (3), the w can be expressed as follows:

w = \frac{F^{\frac{2}{3}}}{{[\frac{16}{9} {(Ē)}^{2} R]}^{\frac{1}{3}}}

(7)

Taking Eq. (7) into Eq. (6), the mathematical expression between the stiffness of the micro-contact pair and the force F can be expressed as follows:

K_{n} = 2 \cdot {(\frac{9}{16})}^{\frac{1}{6}} \cdot Ē^{\frac{2}{3}} \cdot R^{\frac{1}{3}} \cdot F^{\frac{1}{3}}

(8)

Writing the F in Eq. (8) as the form of surface pressure P, that is give F^1/3 multiply the contact area of micro-contact pair and divide by the same value. The contact area of micro-contact pair is πr², in which r² = R ⋅ w. Then Eq. (8) can be expressed as follows:

K_{n} = 2 \cdot {(\frac{9}{16})}^{\frac{1}{6}} \cdot Ē^{\frac{2}{3}} \cdot R^{\frac{1}{3}} \cdot {(π r^{2} \cdot \frac{F}{π r^{2}})}^{\frac{1}{3}} = 2 \cdot {(\frac{9}{16})}^{\frac{1}{6}} \cdot w^{\frac{1}{3}} \cdot Ē^{\frac{2}{3}} \cdot R^{\frac{2}{3}} \cdot π^{\frac{1}{3}} \cdot P^{\frac{1}{3}}

(9)

The size of contact area of the micro-contact pair is much smaller than the joint surface. Therefore, the normal contact stiffness of the micro-contact pair can be considered as the normal stiffness per unit area of the contact at this position. And the mathematical expression between the normal stiffness per unit area and the surface pressure P in the joint can be abbreviated as an exponential form as follows:

K_{n} = A \cdot P^{n}

(10)

In order to obtain the exact relationship between the normal stiffness coefficient of the unit area and the surface pressure in the joint, the parameters A and n in Eq. (10) are obtained by parameter identification. And the normal stiffness distribution of the joint can be obtained though Eq. (10) and the pressure distribution of the joint surface.

b. Relationship between the joint surface pressure and the tangential stiffness per unit area.

The bulge on the joint surface will have tangential deformation when the joint surface which applied the preload force is subjected to the action of the tangential force. According to the literature,¹⁸ the tangential deformation at a single bulge can be expressed as follows

w_{t} = \frac{3 μ F_{n}}{16 \bar{G} r} [1 - {(1 - \frac{F_{t}}{μ F_{n}})}^{\frac{2}{3}}]

(11)

where F_n is the normal load which is caused by the preload force and applied to the micro-contact pair; F_t is the tangential load which is applied to the micro-contact pair; $\bar{G}$ is the equivalent tangential modulus of the micro-contact pair, $1 / \bar{G} = (2 - ν_{1}) / G_{1} + (2 - ν_{2}) / G_{2}$ ⁠; μ is the static friction factor; r is the contact circle radius of the micro-contact pair.

Transforming Eq. (11) into an equivalent form as follows:

F_{t} = μ F_{n} \pm μ F_{n} {(1 - \frac{16 w_{t} \bar{G} r}{3 μ F_{n}})}^{\frac{3}{2}}

(12)

According to Eq. (12), the equivalent tangential stiffness of each micro-contact pair can be expressed as follows:

K_{t} = \frac{d F_{t}}{d w_{t}} = \pm 8 \bar{G} r {(1 - \frac{16 w_{t} \bar{G} r}{3 μ F_{n}})}^{\frac{1}{2}}

(13)

Taking a square on both sides of Eq. (13) to obtain Eq. (14) which can be expressed as follows:

K_{t}^{2} = 64 {\bar{G}}^{2} r^{2} - \frac{1024 {\bar{G}}^{3} r^{3} w_{t}}{3 μ} {(F_{n})}^{- 1}

(14)

Writing the F_n in Eq. (14) as the form of surface pressure P, that is give ${(F_{n})}^{- 1}$ multiply the contact area of micro-contact pair and divide by the same value. The contact area of micro-contact pair is πr², and Eq. (14) can be expressed as follows:

K_{t}^{2} = 64 {\bar{G}}^{2} r^{2} - \frac{1024 {\bar{G}}^{3} r^{3} w_{t}}{3 μ} {(π r^{2} \cdot \frac{F_{n}}{π r^{2}})}^{- 1} = 64 {\bar{G}}^{2} r^{2} - \frac{1024 {\bar{G}}^{3} r w_{t}}{3 μ π} {(P)}^{- 1}

(15)

The size of contact area of the micro-contact pair is much smaller than the joint surface. Therefore, the tangential contact stiffness of the micro-contact pair can be considered as the tangential stiffness per unit area of the contact at this position. And the mathematical expression between the tangential stiffness per unit area and the surface pressure P in the joint can be abbreviated as an exponential form as follows:

K_{t}^{2} = C + B \cdot P^{m}

(16)

where the parameters C, B and m in Eq. (16) are obtained by parameter identification. And the tangential stiffness distribution of the joint can be obtained though Eq. (16) and the pressure distribution of the joint surface.

c. Calculation of stiffness coefficients in Matrix27 stiffness element.

The calculation processes of the stiffness coefficients of all Matrix27 stiffness elements in the bolted joints can be expressed as follows:

Step 1.
In the established finite element model of the assembly specimen structure, applying the bolt preload and making static analysis.
Step 2.
All grid node numbers and its coordinates on the lower block joint surface are selected and represented as set A₁, as well as the data of the upper block is represented as set A₂. Then the pressure values at all grid nodes on lower block joint surface are selected and represented as set A₃.
Step 3.
Classifying the grid node numbers which belong to the same element in the set A₁ and A₂, and expressing as the set B_n, in which n represents the serial number of the set. And set B_n is called as a joint surface element numbered n. The contours of each joint surface element in the joint surfaces of the lower block and the upper block are plotted respectively in MATLAB, which are indicated in FIG. 7(a) and FIG. 7(c). Since the assembly specimen structure is finite element meshed by the Solid187 element, the shape of all joint surface elements on the joint surface is triangular. Compared with the finite element grids at the joint surfaces in ANSYS, which are indicated in FIG. 7(b) and FIG. 7(d), we can find that the grids drawn by MATLAB is in perfect agreement with the grids drawn by ANSYS. And it can be explained that the above classification of grid nodes is correct.
Step 4.
There are three nodes in every joint surface element, so the distribution of the pressure field in joint surface element numbered n can be expressed as follows:

P_{n} (x, y) = P_{00 - n} + P_{01 - n} \cdot x + P_{11 - n} \cdot y

(17)

where x and y are the horizontal coordinate and vertical coordinate at any point in the joint surface element numbered n, P_n(x, y) is the pressure value at any point in the joint surface element numbered n. And based on the coordinates of the grid nodes in each joint surface element and the surface pressure, the parameter P_00-n, P_01-n and P_11-n in the formula are identified.

Step 5.
Based on the grid node coordinates, the mathematical expressions of each contour line in each joint surface element are obtained, as shown in FIG. 8. In which f_n1, f_n2 and f_n3 are the function expressions of the contour lines in the joint surface element which is numbered n. Furthermore, the average pressure of each joint surface element can be obtained by Eq. (18) and Eq. (19).

P_{my - n} = \frac{\int_{y_{2}}^{y_{1}} \int_{\frac{y_{2} - b_{2}}{a_{2}}}^{\frac{y_{3} - b_{3}}{a_{3}}} P_{n} (x, y) d x d y + \int_{y_{1}}^{y_{3}} \int_{\frac{y_{1} - b_{1}}{a_{1}}}^{\frac{y_{3} - b_{3}}{a_{3}}} P_{n} (x, y) d x d y}{S_{n}}

(18)

S_{n} = \int_{y_{2}}^{y_{1}} \int_{\frac{y_{2} - b_{2}}{a_{2}}}^{\frac{y_{3} - b_{3}}{a_{3}}} 1 d x d y + \int_{y_{1}}^{y_{3}} \int_{\frac{y_{1} - b_{1}}{a_{1}}}^{\frac{y_{3} - b_{3}}{a_{3}}} 1 d x d y

(19)

where P_my-n is the average pressure of the joint surface element which is numbered n, P_n(x, y) is the pressure field distribution function of the same joint surface element, S_n is the area of the same joint surface element.

Step 6.
Taking the average surface pressure of each joint surface element into Eq. (10) to figure out the normal stiffness per unit area in the joint surface element, and taking it into Eq. (16) to figure out the tangential stiffness per unit area in the joint surface element. And the total normal stiffness and the total tangential stiffness of the joint surface element can be obtained by multiplying the area of the joint surface element with the normal stiffness per unit area and the tangential stiffness per unit area respectively.
Step 7.
Dividing the total stiffness of a joint surface element equally among all grid nodes on the contour of the joint surface element. Because the same grid node is shared by a number of joint surface elements, the total stiffness of each grid node can be figured out by Eq. (20) which is shown as follows:

K_{M} = \sum_{n = 1}^{N} k_{M - n}

(20)

where K_M is the total stiffness of the grid node which is numbered M, n is the joint surface element number, N is the number of joint surface elements that have the common node M. Thus, the normal stiffness distribution and tangential stiffness distribution based on the joint pressure distribution can be obtained.

Step 8.
Taking the normal stiffness value and the tangential stiffness value of each grid node into the Matrix27 stiffness element at the same position.

FIG. 7.

View large Download slide

Mesh grid distribution map of joint surface. (a) Mesh grid distribution map of upper joint surface. (b) Mesh grid distribution map of upper joint surface in ANSYS. (c) Mesh grid distribution map of lower joint surface. (d) Mesh grid distribution map of lower joint surface in ANSYS.

FIG. 8.

View large Download slide

Oblique edge mathematical expression of element.

The calculation process of stiffness coefficients in Matrix27 stiffness element as shown in FIG. 9.

FIG. 9.

View large Download slide

Flow chart of calculation process of stiffness coefficients in Matrix27 stiffness element.

d. Description of the damping characteristic of the assembly specimen structure.

The damping characteristic of the assembly specimen structure is determined by the structural damping of the components and the damping of the bolted joints. In this paper, the damping characteristic of the assembly specimen structure is described by Rayleigh damping model which is shown as follows:

[C] = α [M] + β [K]

(21)

where [M] is the mass matrix of finite element model of the assembly specimen structure, [K] is the stiffness matrix of finite element model of the assembly specimen structure. And the parameters α and β in the model are obtained by parameter identification.

4. Constraint state

The modal information extracted by the model in free constraint state is more abundant than fixed constraint state, so the constraint state is free constraint in the process of finite element modeling of the assembly structure. On the other hand, the progress of simulating fixed constraint state may generate error when using fixed constraint state. The error will be integrated into the error caused by the dynamic modeling of the joint, which will impact the validity of the proposed joint modeling method. Therefore, the free constraint state is chosen in this research.

II. DYNAMICS EXPERIMENT OF THE BOLTED JOINTS

A. Design of the experiment

As shown in Figs. 10, a dynamic testing device is set up, which includes the PCB’s 086C03 type hammer, two 356A16 type three-way piezoelectric acceleration sensors, a set of LMS’s SCM202 type 16-channel data acquisition instrument, and a set of LMS Test.lab vibration test and analysis software system. The test specimen is formed by the upper block and the lower block, and the assembly specimen structure are connected by four uniform distributed M8 bolts. Besides, the preload torque of each bolt is 25 N.m.

FIG. 10.

View large Download slide

Experimental setup.

In the process of the modal test, the assembly specimen structure is hung on a frame by nylon rope. Because the modal test in this study is to verify the reasonableness and accuracy of the joint dynamic modeling method proposed in this paper. Therefore, the constraint imposed on the assembly structure in the test and the constraint imposed in the dynamic modeling process should be consistent, and the constraint is free constraint. The experimental tools are the PCB’s 086C03 type hammer with nylon hammerhead and 356A16 type three-way piezoelectric acceleration sensor. And the data is collected by the LMS’s SCM202 type 16-channel data acquisition instrument. The first four-order experimental modal frequencies of the assembly specimen is shown in TABLE III, which is extracted by the LMStest.Lab software. The first four-order experimental mode shapes of the assembly specimen structure is shown in TABLE V. The frequency response curves and theirs phase diagrams between all the excitation points and the test points are shown in FIG. 11, which is obtained through the experiment.

FIG. 11.

View large Download slide

Amplitude and phase diagram of frequency response function.

B. Parameter identification of the bolted joint

In this paper, the parameters A, n, C, B, m, α and β in Eq. (10), Eq. (16) and Eq. (21) are identified by the Powell optimization algorithm, which uses the property that the conjugate direction can accelerate the convergence speed. The method can be applied when the derivative of the objective function is not continuous because the it does not need to find the derivative of the objective function. Therefore, the Powell optimization algorithm is a very effective direct search method. This method can be used to solve general unconstrained optimization problems, and can obtain satisfactory results for the objective function optimization problem when the dimension n less than 20.

1. Pressure distribution of the bolted joint

The pressure distribution of the joint surface is the basis for solving the stiffness distribution and parameter identification of the joint surface. In this paper, the static module in ANSYS is used to obtain the pressure distribution of the joint surface. The assembly specimen structure is formed by the upper block and the lower block, and the assembly specimen structure are connected by four uniform distributed M8 bolts. Besides, the preload torque of each bolt is 25 N.m. The preload force for each bolt is

F_{0} = \frac{T}{K d} = \frac{T}{0.2 d} = \frac{25 N \cdot m}{0.2 \times 0.008 m} = 15625 N

(22)

where T is the preload torque of the bolt, d is the nominal diameter of the bolt, F₀ is the bolt preload force. By consulting the mechanical design manual according to the actual installation conditions of the bolts, the value of K is 0.2. Appling the preload force F₀ at the bolt position in the finite element model of the assembly specimen structure, and making static analysis to obtain the pressure distribution of the joint surface, which is shown in FIG. 12.

FIG. 12.

View large Download slide

Pressure distribution on joint surface.

2. Process of the parameter identification

The optimization objective function is established as follows:

f = \sqrt{\sum_{i = 1}^{4} {(ω_{i}^{'} - ω_{i})}^{2}}

(23)

where $ω_{i}^{'}$ is the i-order modal frequency which obtained by theoretical calculations, ω_i is the i-order modal frequency which obtained by experiment, f is the objective function value. When the objective function value is less than the pre-given ε value, the optimization process is finished. The method of determining the value of ε is shown as follows:

ε = \sqrt{\sum_{i = 1}^{4} {(ω_{i} \times 5 %)}^{2}}

(24)

When the Powell optimization algorithm is used to optimize the multi variables, only one variable can be optimized each time. And the process of parameter identification is shown as follows:

Step 1.
Assigning initial values to the unknown parameters (they are the final values of the last round of optimization). And denoting all of the initial values as vector x₀.
Step 2.
Optimizing every variable respectively, and only one variable can be optimized each time.
- The interval of the optimal value of the variable is obtained by the advance and retreat method,²¹ and is expressed by e₁.
- The Fibonacci method is used to reduce the interval by n times to obtain the interval e₂.
- Selecting the median value of the interval as the optimal value of the variable.
Step 3.
Repeating the step 2 until all variables are optimized. And all of the identified parameter values are denoted as vector x_n, where n means it is n round optimization.
Step 4.
Optimizing along the direction from the initial point to the final point of the n round optimization, the optimal values can be obtained and represented by vector y_n.
Step 5.
Taking y_n into Eq. (23), and figuring out the objective function value of the n round optimization. Furthermore, denoting the objective function value as f_n, and Comparing the value of f_n and ε. If f_n < ε, the optimization can be finished. If not, repeating the steps 1-5 until f_n < ε.

Writing an optimize control program in MATLAB which can achieve interactive computing with ANSYS. In this way, MATLAB can automatically call ANSYS, And the initial value of the stiffness of the joint and the initial value of the machine damping coefficient are transmitted to ANSYS for modal analysis. The modal frequency is fed back into the MATLAB cycle control program via the APDL command when the analysis of ANSYS is completed, and the process is automatically iterated until the optimize cut-off condition is reached.

The identified parameters are shown in TABLE II. Taking the parameter values and the joint surface pressure into Eq. (10) and Eq. (16) to figure out the normal stiffness distribution and the tangential stiffness distribution of the joint surface. The pressure distribution, the normal stiffness distribution and the tangential stiffness distribution of the joint in assembly specimen structure are shown in FIG. 13, FIG. 14 and FIG. 15. Furthermore, it can be found that the stiffness near the bolt holes is obviously larger than anywhere else, due to the pressure near the bolt holes is concentrated and it is caused by the preload force of the bolts.

TABLE II.

Value of parameter identification.

A	n	B	C	m	α	β
1.8613e⁵	1.0471	1.5764e¹⁰	2.3786e¹⁰	1.9	0.924	9.2e⁻⁶

FIG. 13.

View large Download slide

Pressure distribution on joint surface.

FIG. 14.

View large Download slide

Normal stiffness distribution on joint surface.

FIG. 15.

View large Download slide

Tangential stiffness distribution on joint surface.

C. Effectiveness verification of the dynamics model of bolted joint

1. Comparison with the experiment

Combining the finite element model of the bolted joint and the finite element model of the upper and lower blocks of the assembly specimen into the integrated finite element model. The modal analysis of the integrated finite element model is carried out by the ANSYS modal analysis, and the calculated results of the first four-order modal frequencies and mode shapes are shown in TABLE III and TABLE V.

TABLE III.

Comparison between theoretical modal values and experimental modal values.

	First order modal(Hz)	Second order modal(Hz)	third order modal(Hz)	Fourth order modal(Hz)
Experimental value	2506.9	3044.5	3746.5	3991.1
The theoretical	2546.4	3173.8	3728.4	3983.9
model in this paper
Error	1.6%	4.2%	0.5%	0.2%

	First order modal(Hz)	Second order modal(Hz)	third order modal(Hz)	Fourth order modal(Hz)
Experimental value	2506.9	3044.5	3746.5	3991.1
The theoretical	2546.4	3173.8	3728.4	3983.9
model in this paper
Error	1.6%	4.2%	0.5%	0.2%

Through comparing the first four-order theoretical modal frequencies and the first four-order experimental modal frequencies of the assembly specimen in TABLE III, it is found that the relative errors are between 0.2% to 4.2%. And the greatest relative error appeared at the second order modal. Furthermore, comparing the first four-order theoretical mode shapes and the first four-order experimental mode shapes of the assembly specimen, which is shown in TABLE V. The comparison results show that the first four-order mode shapes of theoretical model are in basic agreement with those of the experimental method while the two neighboring mode shapes of the theoretical and experimental methods may exchange with each other sometimes. Therefore, the dynamic modeling method which is proposed in this paper is reasonable and feasible.

2. Comparison of results with other researchers

Comparing the results of this paper with the thread-connected spring-damping model in the literature,¹⁹ that is the Yoshimura model in the literature.²⁰ According to the method in literature,²⁰ a spring-damping element is created at each bolt position in the tool holder joint to describe the dynamic characteristics of the bolted joints. And the calculation process of the stiffness coefficient in the spring-damping element is described as follows:

Step 1.
The average surface pressure of the joint surface can be calculated as follows:

P = \frac{4 \times F_{0}}{L \times W} = \frac{4 \times 15625 N}{0.1 m \times 0.14 m} = 4.464285 M p a

(25)

where L is the length of the bolted joints, W is the width of the bolted joints, P is the average surface pressure on the bolted joint surface. And the unit for the dimensions of the joint surface is mm.

Step 2.
Based on the joint surface pressure, the stiffness of the unit area of the joint surface are obtained by the method of the literature.¹⁹

k_{n} = 1.6837 \times e^{11} (N \cdot m) / m^{2}

k_{t} = 8.9728 \times e^{7} (N \cdot m) / m^{2}

where k_n is the normal stiffness per unit area on the joint surface, k_t is the tangential stiffness per unit area on the joint surface.

Step 3.
The integrated normal and tangential stiffness of the bolted joints can be calculated as follows:

K_{n} = k_{n} \times S = 1.6837 \times e^{11} \times 0.014 = 2.35718 \times e^{9} N \cdot m

K_{t} = k_{t} \times S = 8.9728 \times e^{7} \times 0.014 = 1.25619 \times e^{6} N \cdot m

where K_n is the normal total stiffness on the joint surface, K_t is the tangential total stiffness on the joint surface.

Step 4.
Dividing the integrated stiffness of the joint into each spring element in the bolted joints equally, and the stiffness value of each spring element can be calculated as follows:

K_{n}^{'} = K_{n} / 4 = 5.893 \times e^{8} N \cdot m

K_{t}^{'} = K_{t} / 4 = 3.14 \times e^{5} N \cdot m

where $K_{n}^{'}$ is the normal stiffness of a single spring element on the joint surface, K_t is the tangential stiffness of a single spring element on the joint surface The calculated results of the first four-order modal frequencies by the Yoshimura model are shown in TABLE IV, and the first four-order mode shape is shown in TABLE V. Through comparing the first four-order theoretical modals and the first four-order experimental modals of the assembly specimen, it is found that the relative errors are between 26.9% and 43.5%. However, the first four-order theoretical mode shapes obtained by the Yoshimura model are not similar to experimental modal. It is shown that the dynamic modeling method of the joint proposed in this paper is reasonable and feasible, and the joint model that considered the joint pressure model proposed in this paper is more advanced and effective.

TABLE IV.

Comparison between modal values of Yoshimura model and experimental modal values.

	First order modal (Hz)	Second order modal(Hz)	third order modal(Hz)	Fourth order modal(Hz)
Experimental value	2506.9	3044.5	3746.5	3991.1
Yoshimura model	1832.1	1954.6	2117.9	2292.9
Error	26.9%	35.8%	43.5%	42.6%

	First order modal (Hz)	Second order modal(Hz)	third order modal(Hz)	Fourth order modal(Hz)
Experimental value	2506.9	3044.5	3746.5	3991.1
Yoshimura model	1832.1	1954.6	2117.9	2292.9
Error	26.9%	35.8%	43.5%	42.6%

TABLE V.

Comparison between theoretical mode shape and experimental mode shape.

Order	Experimental mode shape	Yoshimura mode shape	The mode shape of the theoretical model in this paper
1
2
3
4

Order	Experimental mode shape	Yoshimura mode shape	The mode shape of the theoretical model in this paper
1
2
3
4

III. CONCLUSION

A new method for dynamic modeling of the bolted joints is proposed in this paper. The uneven pressure distribution on the bolted joint surface influences the modeling accuracy during the dynamic modeling of the bolted joints, and it is considered in this method. Based on the Hertz contact theory, the analytic formulas between the normal and tangential stiffness per unit area and the surface pressure on the joint surface are deduced.
The validity of the proposed dynamic modeling method for the bolted joints is verified by experiments. Comparing the first four-order theoretical modal frequencies by the proposed method and the first four-order experimental modal frequencies, the results show that the relative error is 0.2% to 4.2%. Besides, the first four-order theoretical mode shapes and the first four-order experimental mode shapes are similar and one-to-one correspondence. However, the relative error between the first four-order theoretical modal frequencies and the first four-order experimental modal frequencies is large, which is calculated based on Yoshimura model, and the first four-order theoretical mode shapes and the first four-order experimental mode shapes are not similar.
The dynamic modeling method proposed in this paper can provide a theoretical basis for the accurate dynamic modeling of the bolted joints in machine tools. Based on this modeling method, the dynamic modeling of the bolted joints can be accomplished without experiments. Besides, it also has the advantages of low modeling cost, high modeling efficiency, and suitable for any type of bolted joint.

ACKNOWLEDGMENTS

This research is supported by the National Natural Science Foundation of China under Grant No. 51775452, the National Natural Science Foundation of China under Grant No. 51275426 and the Fundamental Research Funds for the Central Universities under No. A0920502051723-9.

REFERENCES

1.

M.

Burdekin

,

N.

Back

, and

A.

Cowley

,

Mech. Eng. Sci.

21

,

25

(

1979

).

https://doi.org/10.1243/jmes_jour_1979_021_006_02

Google Scholar

Crossref

2.

C. F.

Beards

,

Shock Vib. Dig.

14

,

9

(

1982

).

https://doi.org/10.1177/058310248201400604

Google Scholar

Crossref

3.

G. P.

Zhang

,

Y. M.

Huang

,

W. H.

Shi

, and

W. P.

Fu

,

Int. J. Mach. Tool. Manu.

43

,

699

(

2003

).

https://doi.org/10.1016/s0890-6955(03)00026-9

Google Scholar

Crossref

4.

J. A.

Greenwood

and

J. B. P.

Williamson

,

P. Roy. Soc. Lon. A-Math. Phy.

295

,

300

(

1966

).

https://doi.org/10.1098/rspa.1966.0242

Google Scholar

Crossref

5.

A.

Majumdar

and

B.

Bhushan

,

J. Tribol.

113

,

1

(

1991

).

https://doi.org/10.1115/1.2920588

Google Scholar

Crossref

6.

M.

Yoshimura

,

CIRP Annals

28

,

241

(

1979

).

Google Scholar

7.

S.

Taylor

and

S. A.

Tobias

,

Lumped-constants method for the prediction of the vibration characteristics of machine tool structures

(

Proc. 5th Int. MTDR Conf.

,

1964

).

Google Scholar

8.

C.

Xu

,

J.

Zhang

,

Z.

Wu

,

D.

Yu

, and

P.

Feng

,

Int. J. Adv. Manu.

67

,

1517

(

2013

).

https://doi.org/10.1007/s00170-012-4586-1

Google Scholar

Crossref

9.

H.

Ahmadian

and

H.

Jalali

,

Mech. Syst. Signal Pr.

21

,

1041

(

2007

).

https://doi.org/10.1016/j.ymssp.2005.08.015

Google Scholar

Crossref

10.

L.

Li

,

A. J.

Cai

,

L. G.

Cai

,

X. G.

Ruan

, and

T. N.

Guo

,

J. Vib. Shock

33

,

15

(

2014

).

Google Scholar

11.

J.

Zhang

and

Z. F.

Tong

,

J. Vib. Shock

13

,

15

(

1994

).

Google Scholar

12.

H.

Tian

,

B.

Li

,

H.

Liu

,

K.

Mao

,

F.

Peng

, and

X.

Huang

,

Int. J. Mach. Tool. Manu.

51

,

239

(

2011

).

https://doi.org/10.1016/j.ijmachtools.2010.11.004

Google Scholar

Crossref

13.

K.

Mao

,

B.

Li

,

J.

Wu

, and

X.

Shao

,

Int. J. Mach. Tool. Manu.

50

,

156

(

2010

).

https://doi.org/10.1016/j.ijmachtools.2009.10.017

Google Scholar

Crossref

14.

N. M.

Salih

and

M. J.

Patil

,

Int. J. Adv. Eng. Tech.

3

,

213

(

2012

).

Google Scholar

15.

Y. M.

Huang

,

W. P.

Fu

,

L. X.

Dong

, and

J. X.

Tong

,

Chin. J. Mech. Eng.-En.

29

,

74

(

1993

).

Google Scholar

16.

X. L.

Zhang

,

The Dynamic Characteristics of Machine Joint Surface and Applications

(

Chin. Sci. Tech. Press

,

2002

).

Google Scholar

17.

S.

Wang

and

K. A.

Komvopoulos

,

J. Tribol.

116

,

812

(

1994

).

https://doi.org/10.1115/1.2927338

Google Scholar

Crossref

18.

L.

Kogut

and

I.

Etsion

,

J. Appl. Mech.

69

,

657

(

2002

).

https://doi.org/10.1115/1.1490373

Google Scholar

Crossref

19.

B. Y.

Liao

,

X. M.

Zhou

, and

Z. H.

Yin

,

Modern mechanical dynamics and its application: modeling, analysis, simulation, modification, control, optimization

(

Chin. Mach. Press

,

2004

).

Google Scholar

20.

M.

Yoshimura

,

Mach. Tool.

1

,

142

(

1979

).

Google Scholar

21.

C. F.

Ma

,

Optimization method and its Matlab program design

(

Sci. Press

,

Beijing

,

2010

).

Google Scholar

2018

Author(s)

Dynamic modeling method of the bolted joint with uneven distribution of joint surface pressure

I. A new dynamic model of the bolted joints in machine tools

A. Description of the research object

B. Dynamic modeling of the bolted joints

1. Mesh generation

2. Establishing the Matrix27 stiffness elements on the bolted joints

3. Calculation of stiffness and damping coefficient

a. Relationship between the joint surface pressure and the normal stiffness per unit area.

b. Relationship between the joint surface pressure and the tangential stiffness per unit area.

c. Calculation of stiffness coefficients in Matrix27 stiffness element.

d. Description of the damping characteristic of the assembly specimen structure.

4. Constraint state

II. DYNAMICS EXPERIMENT OF THE BOLTED JOINTS

A. Design of the experiment

B. Parameter identification of the bolted joint

1. Pressure distribution of the bolted joint

2. Process of the parameter identification

C. Effectiveness verification of the dynamics model of bolted joint

1. Comparison with the experiment

2. Comparison of results with other researchers

III. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Dynamic modeling method of the bolted joint with uneven distribution of joint surface pressure

I. A new dynamic model of the bolted joints in machine tools

A. Description of the research object

B. Dynamic modeling of the bolted joints

1. Mesh generation

2. Establishing the Matrix27 stiffness elements on the bolted joints

3. Calculation of stiffness and damping coefficient

a. Relationship between the joint surface pressure and the normal stiffness per unit area.

b. Relationship between the joint surface pressure and the tangential stiffness per unit area.

c. Calculation of stiffness coefficients in Matrix27 stiffness element.

d. Description of the damping characteristic of the assembly specimen structure.

4. Constraint state

II. DYNAMICS EXPERIMENT OF THE BOLTED JOINTS

A. Design of the experiment

B. Parameter identification of the bolted joint

1. Pressure distribution of the bolted joint

2. Process of the parameter identification

C. Effectiveness verification of the dynamics model of bolted joint

1. Comparison with the experiment

2. Comparison of results with other researchers

III. CONCLUSION

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Submit your article

Sign up for alerts

Related Content

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only