Nonlinear Vibration Induced by Friction in a Ball Joint System

Abstract

In this paper, the nonlinear oscillations induced by friction in a ball-on-socket system are investigated. The nonlinear time response was obtained by solving the differential equations of the friction-noise model of the finite element ball with multiple modes. The different patterns of motion were analyzed via the bifurcation diagram, Poincare map, and recurrence plot. The Lyapunov exponents of the discontinuous system with distributed contact were calculated using the Muller method. From the analysis, it is shown that a friction-noise of a ball joint can retain periodic, quasi-periodic, or chaotic oscillations with respect to tilted contact.

1. Introduction

A ball joint is a major mechanical part consisting of a spherical bearing and a semi-spherical socket. It is widely used in the suspension and steering linkages of automotive, robotic arms, and even biomechanics such as human hip joints. Fretting and squeaking is one frequent issue; thus, consistent lubrication is essential for the prevention of these problems. The friction-induced vibration is the self-excited vibration [1,2]. Its mechanisms have long been investigated in the application of automotive brake and gear systems. The linear stability analysis has been adopted for estimating the onset of squeal [3,4,5,6,7,8].

The numerical simulation for friction noise in a ball-on-socket has been also carried out on the basis of complex eigenvalue analysis [9,10,11,12]. It has been suggested that the positive real parts of the system’s eigenvalues indicate the unstable modes generating friction noise frequencies. However, the results of the linear stability theory did not properly estimate the squeak frequencies in the experiment. Therefore, the transient time-domain analysis should be considered to reproduce the experimental squeak frequencies, where the nonlinear differential equations subject to friction sliding are numerically solved. Recently, Kang’s nonlinear finite element (FE) model qualitatively reproduced the frequency spectrum in an experiment [13].

It is known that the friction-induced vibrating system has sufficient complexity. It may contain chaotic motions as well as regular ones. Simple analytical models have been proposed to investigate the rich pattern of nonlinear oscillations [14,15,16]. The Lyapunov exponents can indicate the sensitivity to the initial values as the characteristics of chaos in the manner that the nonlinear motion exhibits the chaotic behavior when the largest Lyapunov exponent is positive. In real applications, the degree-of-freedom in the FE model is tremendously increased. In this case, the time-delay reconstruction for the time series data of the experiment or simulation was carried out to calculate the Lyapunov exponents [17,18,19]. The recurrence quantification analysis (RQA, [20]) was also used to display the recurrence pattern of the time series data in brake squeal application [21,22]. Recently, CNN-based deep learning was implemented to determine the classification of recurrence patterns of friction-induced vibrations [23,24].

To the author’s knowledge, this is the first work to present the calculation of the spectrum of Lyapunov exponents for the distributed FE friction contact on the ball joint system in order to determine the chaotic behavior of ball joint squeak. Then, it clarifies that the frequency spectrum can be quite variable to small changes in the parameters such as the tilt angle of axial load. This approach can be applied to hip joint squeak problems with a complex geometry and various frequency spectrums in biomedical applications.

We also attempt to seek the different types of attractors for the nonlinear oscillation of the FE ball-on-socket system in the analytical and numerical manner. In particular, the possibility of a chaotic attractor in the ball joint system is proposed through nonlinear analysis such as a bifurcation diagram, Poincare map and recurrence analysis. The degree of freedom in the FE ball joint model was reduced through the truncated Galerkin approximations to the original nonlinear equations of motion. Then, the Lyapunov exponents were directly estimated by solving the differential equations. We also clarify that the unstable modes in the linear stability analysis do not predict the squeak frequencies of the transient time-domain analysis under several circumstances.

2. Materials and Methods

The ball joint system was constructed in the finite element manner as shown in Figure 1a. The ball joint rotated in the $x$, $y$, and $z$ axes when the tilted reaction pre-load was applied on the center of the ball. The tilting angle between the axis of the beam and reaction pre-load was considered to be the system parameter. The clearance between the ball and socket produced the contact area and the contact pre-stress distribution as in Figure 1b.

The highly nonlinear contact kinematics were formulated on the contact between the ball and socket. Contact area on the ball was estimated by Hertz contact theory [10,11,12,13]. The analytical contact area was projected onto the FE geometry; thus, the node within the FE contact surface was called the contact node. In the perturbed state from the steady-sliding equilibrium, the contact force vector normal to the contact surface was defined on the $k$th contact node as:

$$N_{\mathrm{c}}^k=-N_{\mathrm{c}}^k{e}_r$$

(1)

where the normal force is the summation of the normal pre-stress $p_{\mathrm{n}}^k$ and the normal load variation $k_{\mathrm{c}}{u}_r^k$ if the $k$th contact node is under contact, and zero if it is non-contacted. So, it can take the following the form:

$$N_{\mathrm{c}}^k=\left\{p_{\mathrm{n}}^k+k_{\mathrm{c}}\left(u_r^k-u_{r,\mathrm{eq}}^k\right)\right\}\cdot \Theta \left\{p_{\mathrm{n}}^k+k_{\mathrm{c}}\left(u_r^k-u_{r,\mathrm{eq}}^k\right)\right\}$$

(2)

where the subscript ‘eq’ denotes the equilibrium state and $\Theta (\cdot )$ is the Heaviside function representing the discontinuous property of contact modelling such that $\Theta (x)=0$ if $x<0$, and $\Theta (x)=1$ otherwise.

The friction force at the $k$th contact node is described in the Coulomb’s friction law as:

$$f_{\mathrm{c}}^k=-{\mu }_{}^k{N}_{\mathrm{c}}^k\frac{V_{\mathrm{c}}^k}{\left|V_{\mathrm{c}}^k\right|}$$

(3)

where the velocity vector of the $k$th contact node is the time-derivative of the perturbed position vector as:

$$V_{\mathrm{c}}^k=\frac{\mathrm{d}}{\mathrm{d}t}\left(x_{\mathrm{c}}^k+u_{\mathrm{c}}^k\right)$$

(4)

where $x_{\mathrm{c}}$ is the position vector from the ball center and $u_{\mathrm{c}}$ is the deflection vector of the contact node. The friction coefficient is described in the smoothing curve method [25,26] as:

$${\mu }_{}^k=\left(1-e^{-\mathrm{d}\left|V_{\mathrm{c}}^k\right|}\right)\left\{{\mu }_{\mathrm{d}}+\left({\mu }_{\mathrm{s}}-{\mu }_{\mathrm{d}}\right)\cdot e^{-\mathrm{f}\left|V_{\mathrm{c}}^k\right|}\right\}$$

(5)

where $\mathrm{d}$ and $\mathrm{f}$ are the slope factors in the creep and sliding regions of the friction curve, respectively, and ${\mu }_{\mathrm{d}}$ and ${\mu }_{\mathrm{s}}$ are the kinetic and static friction coefficients of the friction curve.

The virtual work performed by the contact normal and friction forces is then defined and numerically estimated in the finite element manner as:

$$\delta {\mathrm{W}}_{\mathrm{c}}={\displaystyle {\int }_{{\mathrm{A}}_{\mathrm{c}}}\left(f_{\mathrm{c}}+N_{\mathrm{c}}\right)}·\delta u_{\mathrm{c}}\mathrm{dA}\cong {\displaystyle \sum _{k=1}^{{\mathrm{N}}_{\mathrm{nd}}}\left(f_{\mathrm{c}}^k+N_{\mathrm{c}}^k\right)·\delta u_{\mathrm{c}}^k\Delta {\mathrm{A}}^k}$$

(6)

where ‘$·$’ denotes the vector inner product, ${\mathrm{N}}_{\mathrm{nd}}$ is the number of contact nodes, $\Delta {\mathrm{A}}^k$ is the finite area of the $k$th contact element, and the contact deflection vector $u_{\mathrm{c}}$ is expressed in the modal expansion form:

$$u\left(x_{\mathrm{c}}^k,t\right)={\displaystyle \sum _{i=1}^{\mathrm{N}}\left[{\Phi }_{\theta ,i}^k{e}_{\theta }+{\Phi }_{\varphi ,i}^k{e}_{\varphi }+{\Phi }_{r,i}^k{e}_r\right]}\cdot q_i(t)$$

(7)

where $\mathrm{N}$ represents the truncated number of the vibration modes, $q_i$ is the $i$th modal coordinate, and ${\Phi }_{\theta ,i}$, ${\Phi }_{\varphi ,i}$, and ${\Phi }_{r,i}$ are the elements of the $i$th eigenvector in the $\theta$, $\varphi$, and $r$ directions of the spherical coordinates, respectively.

The discretized nonlinear equations of motion with distributed and discontinuous FE friction contact can reduce to:

$${\ddot{q}}_i+{\omega }_i^2{q}_i={\displaystyle \sum _{k=1}^{{\mathrm{N}}_{\mathrm{nd}}}\left\{\begin{array}{c}-{\mu }^k\frac{N_{\mathrm{c}}^k}{\left|V_{\mathrm{c}}^k\right|}\left(V_{\mathrm{c}}^k·e_{\theta }\right){\Phi }_{\theta ,i}^k\\ -{\mu }^k\frac{N_{\mathrm{c}}^k}{\left|V_{\mathrm{c}}^k\right|}\left(V_{\mathrm{c}}^k·e_{\varphi }\right){\Phi }_{\varphi ,i}^k\\ -{\mu }^k\frac{N_{\mathrm{c}}^k}{\left|V_{\mathrm{c}}^k\right|}\left(V_{\mathrm{c}}^k·e_r\right){\Phi }_{r,i}^k\\ -N_{\mathrm{c}}^k{\Phi }_{r,i}^k\end{array}\right\}}\Delta {\mathrm{A}}^k$$

(8)

where ${\omega }_i^{}$ is the circular natural frequency of the $i$th mode. The linear stability analysis can be conducted by linearizing Equation (8) at the equilibrium.

The Lyapunov spectrum is determined by solving the derived differential equations of motion and the variational equation simultaneously. The nonlinear differential Equation (8) is rewritten in the state-space form as:

$$\dot{z}=F(z),\begin{array}{c}\end{array}z(t_0)=z_0$$

(9)

where $F$ is the vector of the nonlinear functions, and

$$\begin{array}{l}z=\left\{\begin{array}{ccccccc}{z}_1& z_2& \cdots & z_{2k-1}& z_{2k}& \cdots & z_{2\mathrm{N}}\end{array}\right\}\\ =\left\{\begin{array}{ccccccc}{q}_1& {\dot{q}}_1& \cdots & q_k& {\dot{q}}_k& \cdots & {\dot{q}}_{\mathrm{N}}\end{array}\right\}\end{array}$$

(10)

It should be noted that Equation (9) represents the discontinuous dynamical system due to the non-smooth normal stress function in Equation (2) over the distributed contact area. Therefore, the algorithm of calculating the Lyapunov exponents of the non-smooth dynamical system is required.

Once the time and locations at the instance ($t_i$) between the contact and non-contact at every FE node are detected, Equation (9) can be discretized in the time step by Muller’s method [27,28]:

$$t_{i-1}<t<t_i:\dot{z}=F_i(z),\begin{array}{c}\end{array}z(t_{i-1})=z(t_{{}_{i-1}}^{+})$$

(11)

$$t=t_i:0=h(z(t_i^{-}))$$

(12)

$$z(t_i^{+})=g(z(t_i^{-}))$$

(13)

$$t_i<t<t_{i+1}:\dot{z}=F_{i+1}(z),\begin{array}{c}\end{array}z(t_i)=z(t_i^{+})$$

(14)

It is noted that the friction contact is discontinuous because each contact node can experience the transition between the contact and non-contact condition. This transition of each contact node on FE contact model should be carefully dealt with alongside the transition time as:

$$\delta t_i={\overline{t}}_i-t_i=-\frac{{\frac{\partial h}{\partial z}|}_{z=z_i^{-}}\cdot \delta z_i^{-}}{{\frac{\partial h}{\partial z}|}_{z=z_i^{-}}\cdot F_i(z_i^{-})}$$

(15)

$$\delta z_i^{+}={\frac{\partial g}{\partial z}|}_{z=z_i^{-}}\cdot \delta z_i^{-}+\left[{\frac{\partial g}{\partial z}|}_{z=z_i^{-}}\cdot F_i(z_i^{-})-F_{i+1}(z_i^{+})\right]\delta t_i$$

(16)

where $z_i^{-}=z(t_i^{-})$ and $z_i^{+}=z(t_i^{+})$, the plus and minus represent the right- and the left-sided limits, $h$ and $g$ are the indicator function and transition condition, respectively, and ${\frac{\partial g}{\partial z}|}_{z=z_i^{-}}=I$ for the piecewise linear stiffness model. Moreover, the indicator function $h$ is discretized for all FE contact nodes in Equation (2) and the discontinuity condition in the distributed contact geometry can be detected by applying the single equation form representing the indicator of transition over all contact nodes between the ball and socket as follows:

$$h(t_i)={\displaystyle \prod _{j=1}^{n_{\mathrm{c}}}{h}_j(t_i)}=s_1(t_i)\cdot s_2(t_i)\cdots s_{n_{\mathrm{c}}}(t_i)=0,$$

(17)

$$h_j(t_i)=\left\{p_{\mathrm{n}}^k+k_{\mathrm{c}}\left(u_r^k-u_{r,\mathrm{eq}}^k\right)\right\}\cdot \Theta \left\{p_{\mathrm{n}}^k+k_{\mathrm{c}}\left(u_r^k-u_{r,\mathrm{eq}}^k\right)\right\}=s_j(t_i),\begin{array}{c}\end{array}j=1,2,\cdots ,n_{\mathrm{c}}$$

(18)

The corresponding Jacobian matrix of the indicator function is numerically estimated as:

$$\frac{\partial h}{\partial z}=\left\{\begin{array}{cccccccc}\frac{\partial h}{\partial z_1}& \frac{\partial h}{\partial z_2}& \cdots & \frac{\partial h}{\partial z_{2k-1}}& \frac{\partial h}{\partial z_{2k}}& \cdots & \frac{\partial h}{\partial z_{2\mathrm{N}}}& \frac{\partial h}{\partial z_{2\mathrm{N}+1}}\end{array}\right\}$$

(19)

where $\partial h/\partial z_{2n}=0,\begin{array}{c}\end{array}n=1,2,\cdots ,\mathrm{N}$.

The disturbed trajectory is described in the form:

$$\delta \dot{z}=\frac{\partial F(z(t))}{\partial z}\cdot \delta z,\begin{array}{c}\end{array}\delta z(t_0)=\delta z_0$$

(20)

where the solution of Equation (20) takes the following form:

$$\delta z(t)={\Phi }_t(z_0)\delta z_0$$

(21)

where ${\Phi }_t(z_0)$ is the state transition matrix that is obtained by solving the following variational equation:

$${\dot{\Phi }}_t(z_0)=\frac{\partial F(z(t))}{\partial z}\cdot {\Phi }_t(z_0),\begin{array}{c}\end{array}{\Phi }_{t_0}(z_0)=I$$

(22)

The Lyapunov exponents are defined as:

$${\lambda }_i=\underset{\tau \to \infty }{\lim }\frac{1}{{\tau }_{\mathrm{o}}}\ln \left|m_i({\tau }_{\mathrm{o}})\right|$$

(23)

where $m_1(t)$, $m_2(t)$, …, $m_{2\mathrm{N}}(t)$ are the eigenvalues of ${\Phi }_t(z_0)$ for any initial point $z_0$, and ${\tau }_{\mathrm{o}}=t/{\omega }_1$. Here, the Lyapunov exponents can be numerically calculated by the Gram–Schmidt orthonormalization [18,19], where the largest Lyapunov exponent indicates the sensitivity to the initial values determining the chaotic attractor.

Alternatively, the recurrence analysis can also be conducted for the nonlinear time-series data of the numerical solution. For this, the motion of the central contact point in the sliding direction of the ball joint is numerically calculated by the modal solutions of Equation (9) and the fully-developed time-series are recorded by sampling time $\Delta t$ as:

$$u_s(n\Delta t)=u\left(x_{\mathrm{ctr}},n\Delta t\right)·V_{\mathrm{o}}^k/\left|V_{\mathrm{o}}^k\right|$$

(24)

A time-delay reconstruction in $m$ dimension is then taken for the selected time-series as [29]:

$$y_n=\left(u_{s,n-(m-1)\tau },u_{s,n-(m-2)\tau },\cdots ,u_{s,n-\tau },u_{s,n}\right)$$

(25)

where the embedding dimension $m$ and the time lag $\tau$ are determined by the mutual information method [30] and the false nearest neighbors (FNN) [31]. For the recurrence analysis [20], the recurrence plot (RP) is firstly defined to measure the recurrences of a trajectory in the reconstructed phase space as:

$$R_{i,j}(\epsilon )=\Theta (\epsilon -\Vert y_i-y_j\Vert )$$

(26)

where $\Vert \cdot \Vert$ is the Euclidean norm and $\epsilon$ is the threshold distance.

3. Results

In this study, the steel beam had a 20 mm diameter and 200 mm length and the steel ball with 31.75 mm diameter was in contact with the aluminum socket with a 31.95 mm diameter. For the loading and movement conditions, the reaction pre-load was 150 N, and the rotation components (${\Omega }_1$, ${\Omega }_2$, ${\Omega }_3$) were ($1/\sqrt{2}$, $1/\sqrt{2}$, $0$) [rad/s]. The eleven vibration modes of the ball joint up to 13 kHz were used in the analysis as shown in Figure 2. In the numerical simulation, the Stribeck-type friction curve under mixed lubrication was applied to the nonlinear ball joint model as illustrated in Figure 3.

It has been shown that the tilting angle of the reaction pre-load force is an important parameter in ball joint squeaking [13]. Thus, the tilting angle was chosen to be the main parameter in this study. The different patterns of motion in ball joint squeaking were demonstrated for the several tilting angles, 2.0°, 3.5°, and 5.0°. The mismatch between the linear stability analysis and the frequency spectrum of the nonlinear time series will also be well illustrated.

For the 2.0° tilting angle, the four unstable modes were detected in the linear stability analysis by determining the complex eigenvalues of the linearized equations of Equation (8) as shown in Figure 4a. The corresponding nonlinear time series and its phase plot are demonstrated in Figure 4b,c for the limit cycle oscillations of the central contact node in the sliding direction. One fundamental frequency of 4.1 k Hz and its higher harmonics in the frequency spectrum of the time series are shown in Figure 4d. It was found that the fundamental frequency corresponded to the second bending mode. Even for the unstable multiple modes, the second bending mode dominated the nonlinear oscillations. The topologically similar phase trajectories were reconstructed as shown in the three-dimensional plot of Figure 4e. Its recurrence pattern can be visualized via the recurrence plot in Figure 4f by applying Equation (26) to the reconstructed time-delay vectors, where the diagonal lines were repeated with an equally spaced period of the fundamental frequency as the periodic attractor.

A similar analysis was conducted for the tilting angle of 3.0° as illustrated in Figure 5. Figure 5a shows that the four unstable modes in the linear stability analysis were similar to those in the case of 2.0°. However, the dominant frequency in the frequency spectrum (Figure 5d) turned out to be 12.48 kHz, which was the fourth bending mode. This implies that the dominant mode in the nonlinear oscillations shifted from the second bending mode to the fourth bending mode at this configuration. The nonlinear time series, its phase trajectories, and the time-delay phase trajectories indicate that this motion was nearly periodic (Figure 5b,c,e). In the enlarged RP of Figure 4f, the diagonal lines were mainly repeated with a period of 1/12.48 kHz.

For the tilting angle of 3.5°, the five unstable modes were found in the linear stability analysis as in Figure 6a. Its oscillation pattern measured at the central contact node of the ball had periodic waviness as in Figure 6b, and the corresponding phase trajectories were not purely periodic (Figure 6c). In the spectrogram of Figure 6d, the frequency of 0.57 kHz ($={\mathrm{f}}_2-3\times {\mathrm{f}}_1$) newly appeared by the superposition of the 3.87 kHz ($\equiv {\mathrm{f}}_1$, the second bending mode) and 12.18 kHz ($\equiv {\mathrm{f}}_2$, the fourth bending mode). Then, two-periodic trajectories were coupled with the 0.57 kHz and 3.87 kHz frequencies. The ratio of the two frequency components was irrational, so the composite motion took a torus as the quasi-period oscillations. This is demonstrated in the three-dimensional time-delay trajectories in Figure 6e. In the RP of Figure 6f, the diagonal lines were mainly repeated with an equally spaced period of 1/0.57 kHz, which is equivalent to the period of the waviness in Figure 6b.

On the other hand, the strange oscillation pattern was observed for the tilting angle of 5.0°. The linear stability analysis resulted in the five unstable modes similar to the above (Figure 7a). However, the oscillations seemed aperiodic as shown in the measured time series (Figure 7b) and its phase trajectories (Figure 7c). The frequency spectrum of this chaotic oscillation had peaks similar to the quasi-periodic case, but it contained the noise-like broadband spikes (Figure 7d). The time-delay trajectories in Figure 7e imply a strange attractor where the corresponding RP did not have any recurrence pattern (Figure 7f).

For the recurrence analysis and Poincare mapping, the time-delay reconstruction was taken for the nonlinear time series by selecting the time delay and embedding dimension from the mutual information and the false nearest neighbors (FNN). The first minimum for the chaotic case of the tilting angle 5.0° in Figure 8a was nearly 20, which was chosen to be the time delay in this analysis. By checking the false nearest neighbors in Figure 8b, the optimal embedding dimension was chosen to be 5; thus, the reconstructed phase trajectories in Figure 4e, Figure 5e, Figure 6e, and Figure 7e were topologically similar to the physical phase plot as in Figure 4c, Figure 5c, Figure 6c, and Figure 7c.

Then, the transition from periodic to quasi-periodic and chaotic oscillations was investigated by the recurrence analysis. Figure 9 illustrates that the successive RPs from 2.4° to 4.6° by the increment of 0.2° were captured within 5 milliseconds. For 2.4°~2.8°, the periodic line structure covered the entire RP nearly with a period of ¼ kHz. As the tilting angle increased up to 3.4°, the period of the periodic lines turned to the reciprocal of the lower frequency component in quasi-periodic oscillation as 1/0.57 kHz. Then, the period of the periodic lines became double for 3.6°~4.0°. Above the value of 4.2°,, the recurrence pattern was not seen any more in this local RP.

The Poincare map corresponding to the successive tilting angles was also illustrated by choosing the cross-section of the time-delay phase trajectories as shown in Figure 10. For the periodic trajectories, the set of narrow points were plotted on the map. For the quasi-periodic trajectories, a single embedded circle was seen in the manner that the cross-section intersected the torus geometry. Then, the set of points on the Poincare section scattered arbitrarily as they behaved as the chaotic attractor.

For a given friction curve in Figure 3 and the differential Equation (8), the Lyapunov exponents in Equation (23) were estimated by directly solving Equation (20) with respect to the tilting angle variation. Figure 11a shows that the largest Lyapunov exponent was nearly zero from 1.0° to 4.0°, but it clearly became positive above the value of 4.0°. Particularly, the second and third Lyapunov exponents were also seen to be positive nearly from 5.0° to 8.0°, which indicates the hyper-chaos. The bifurcation diagram in Figure 11b confirms that the chaotic behavior occurred at the higher tilting angles. It reveals that the nonlinear oscillations of the ball joint may have generated the chaotic signal under the relevant configuration.

4. Conclusions

This paper dealt with the nonlinear motions induced by friction in the ball joint system. The different pattern of the motion was obtained from the nonlinear differential equations based on the FE geometry and analytical contact kinematics. Simultaneously, the unstable modes were calculated from the linearized equations of motion at the equilibrium of the friction-engaged ball joint system. The results are summarized as follows:

• The proposed approach can solve the Lyapunov exponents of the complex FE model with distributed frictional contacts;
• The unstable frequencies of the linear ball joint model were found to mismatch with the ball joint squeak frequencies of the numerical nonlinear response;
• In the nonlinear analysis, one fundamental frequency appeared in the periodic motion, and two incommensurate fundamental frequencies were seen in the quasi-periodic and chaotic motions;
• The tilting angle of the reaction pre-load was shown to be the influential factor determining the oscillation pattern. As the tilting angle increased, the periodic motion turned to quasi-periodic, and eventually, chaotic behavior. It revealed that the frequency spectrum and oscillation of ball joint squeaking can drastically change with the tilting angle.