Chaotic characteristics and attractor evolution of friction noise during friction process

Friction experiments are conducted on a ring-on-disk tribometer, and friction noise produced during the friction process is extracted by a microphone. The phase trajectory and chaotic parameters of friction noise are obtained by phase-space reconstruction, and its attractor evolution is analyzed. The results indicate that the friction noise is chaotic because the largest Lyapunov exponent is positive. The phase trajectory of the friction noise follows a “convergence-stability-divergence” pattern during the friction process. The friction noise attractor begins forming in the running-in process, and the correlation dimension D increases gradually. In the stable process, the attractor remains steady, and D is stable. In the last step of the process, the attractor gradually disappears, and D decreases. The friction noise attractor is a chaotic attractor. Knowledge of the dynamic evolution of this attractor can help identify wear state changes from the running-in process to the steady and increasing friction processes.


Introduction
Relative motion between two bodies causes friction and contact surface wear, which may exert negative effects on the reliability, security and usage of equipment [1,2]. Distinguishing and predicting wear states is thus important to extend the life of a mechanical system. However, identifying the wear states of complicated nonlinear dynamic processes is difficult.
In recent years, many domestic and foreign researchers have published studies on the mechanisms and factors affecting friction noise [3]. Chen et al. [4], for example, performed measurements of the noise induced by friction, and used the waveform and power spectrum to analyze the friction noise sound pressure level. The results of this indicated that friction-induced noise was generated by relative sliding friction and vibration motion. Chen and Zhou [5] applied the concept of friction coefficient and observations of scar topography to analyze the mechanism of friction noise. The results of this work showed that fluctuations in friction force were the main mechanism of this noise. Lars and Staffan [6] studied a spiral-shaped modification of the surface topography of a brake disk to reduce noise; results suggested that a spiral pattern could strongly reduce squealing. A basic study of friction noise caused by fretting was conducted by Jibiki et al. [7], and the results of this research indicated that friction noise was generated by certain cycles of fretting. The sound pressure level increased with increasing fretting stroke and frequency, which is related to the average sliding velocity and the wear loss. Chen et al. [8] studied the relationship between friction noise and friction surface characteristics under reciprocating sliding. The group found that the worn morphology with high intensity level noise presented obvious flake fractures, and that these fractures were important causes of friction noise. Le Bot et al. [9] carried out an experiment on friction noise during dry contact under light pressure, and results demonstrated that friction noise increased with sliding speed and roughness increasing. Disk surfaces with different groove structures were investigated on a pad-disk tester by Wang et al. [10,11], who demonstrated that a disk surface with 45 degree grooves could strongly reduce friction noise.
The main purpose of studying friction signals is to identify wear states in the friction process. However, the friction process is extremely complex, and accurately identifying and predicting wear conditions from previous studies on friction noise is difficult. Therefore, several researchers have sought to determine the nonlinear dynamics of friction signals to recognize wear states. Zhu et al. [12] extracted the friction force and vibration signals from a pin-disk experiment to realize the chaotic characteristics of a tribological system, and utilized the spectrum and fractal dimension methods to study these signals quantitatively. The results of this research showed that friction signals presented chaotic natures, and that the tribological system was a chaotic system. Zhou et al. [13] studied the chaotic characteristics of friction temperature in the friction process and found that the temperature signal featured chaotic characteristics that could be utilized to recognize changes in wear states. Sun et al. [14,15] applied chaos theory to study the friction vibration signals extracted from a pin-disk test, and results revealed that these signals demonstrated chaotic characteristics. Wear states can be identified by the evolution of friction vibration attractors. Liu et al. [16] conducted spherical-on-disk running-in tests, and chaotic attractors were used to analyze cross correlations between tangential and normal vibrations. The results of this work demonstrated that the chaotic attractors of vibration signals converged as running-in continued and could be utilized to describe changes in wear states. Oberst and Lai [17] utilized a recurrence plot to reveal the chaotic characteristics of brake squeal noise.
Friction noise contains information reflecting wear states, and friction noise signals can be collected in real time without affecting the normal friction process. Therefore, nonlinear dynamics theory may be utilized to prove that friction noise presents chaotic characteristics and that chaotic attractors undergo dynamic evolution. Friction noise thus presents a new route through which wear states can be monitored on-line, identified and predicted. Based on the rotational movement of a ring-disk under oil lubrication, the aim of this paper is to illustrate the complex chaotic characteristics and attractor evolution of friction noise signals, knowledge of which is instructive in revealing wear states. In Section 2, the experimental apparatus and specimens are introduced, the tests are conducted, and the original time series of the friction noise signals are obtained. In Section 3, the processed time series of the friction noise signals are obtained by means of reducing sampling and de-nosing. In Section 4, the evolution of the phase trajectories of friction signals is presented by reconstructing the phase space. In Section 5, the correlation dimension D and largest Lyapunov exponent are calculated. Finally, in Section 6, a discussion and analysis are provided, and the main conclusions of this work are given.

Tribometer description
The friction experiments were carried out on a rotating tribometer. Friction noise signals can be extracted in real-time via the ring-disk rotating process. The experimental device is composed of a power system, a loading system, a clamp system, and a data acquisition and an analysis system. The equipment is shown in Fig. 1. The upper sample ring was installed on a ring holder by a locking pin and a locking bolt, while the lower sample disk was mounted on a disk holder by a locating pin illustrated in Fig. 2. The friction torque was measured by a torque sensor, which was also attached to the disk holder. Friction noise signals were measured by a microphone PCB fixed 35 mm from the center of the friction pairs. The sensitivity of the microphone PCB was 45 mV/Pa, and its dynamic range was 15-122 dB. Signals were extracted and stored by a CoCo80 analyzer.

Test samples and test conditions
The test samples included ring and disk specimens. The upper sample was a ring made of GCr15 bearing steel and with a Rockwell hardness of 61 HRC after quenching treatment. Its outside diameter was 34 mm, and its inner diameter was 24 mm. The surface of this ring was first processed by turning, and then ground and polished using sandpapers of 800#, 1200#, 1500#, and 2000# in sequence. The roughness Ra of this specimen was measured, and the mean of three measurements was considered the final Ra of the specimen. The ring specimens showed a roughness of Ra=0.040−0.043 μm. The lower sample was a disk made of AISI 1045; its Rockwell hardness was 44 HRC without heat treatment and its diameter was 46 mm. The disk had an initial roughness of Ra=5.520− 6.000 μm after turning. The equivalent radius of the friction pair was 14.5 mm, and the nominal contact area was 455.53 mm 2 . Table 1 shows the experimental conditions of five tests.
The tests were conducted under atmospheric conditions (at 16−22 °C ; relative humidity, 48%−58%). Prior to testing, the specimens were cleaned with ethanol (97% pure) using an ultrasonic cleaner. A volume of 0.2 ml of 15W-40 lubricating oil was dropped onto the working surface of the lower sample, and good contact of the specimen surfaces was ensured    | https://mc03.manuscriptcentral.com/friction during installation. The sampling frequencies of the sound level sensor and the friction torque were 12.8 kHz and 300 Hz, respectively, and the background noise was measured before the tests. To ensure the repeatability of the test results, five tests were carried out at least three times to account for the randomness of the friction noise.

Experimental signal processing
The voltage signals in the experiments were converted to sound pressure levels in decibels using the equation where, U is the voltage of the test data, U 0 is the voltage of the background noise, P is the sound pressure level of the test data after conversion, P 0 is the sound pressure level of the background noise after conversion, and sty is the sensitivity of the microphone.
where P c is the reference value of sound pressure, P c = 2 × 10 −5 Pa, Lp is the sound pressure level of the test data, and Lp 0 is the sound pressure level of the background noise. The equation used to calculate the sound pressure level after removal of the background noise is Lp 10 log (abs(10 10 )) where Lp 1 is the sound pressure level of the friction noise without the background noise. Generally speaking, the friction noise of the friction system includes other insignificant noise that may influence the dynamic characteristics of the former. Therefore, filtering and denoising of the friction noise signals was conducted by empirical mode decomposition (EMD).
EMD was put forward by Huang et al. [18] to filter nonlinear and unstable time series. Using EMD, time series can be decomposed into a set of finite and intrinsic mode functions (imf) with different components. Then, specific imf components selected from the original components are reset to a new time series to analyze and calculate the data during follow-up work [19,20]. Details of the EMD algorithm are as provided as follows.
The maximum and minimum values of a time series x(t) are determined, and the upper envelope curve x max (t) and lower envelope curve x min (t) of the original time series are calculated by the cubic spline functions. The mean of the upper and lower envelope curves is denoted as m(t). The original time series x(t) minus the mean m(t) yields a new series h(t) without lowfrequency components. The relevant functions are shown in Eqs. (4) and (5).
In general, the first h(t) is not necessarily a stationary data series. Therefore, it should be calculated using Eqs. (4) and (5) repeatedly until it meets the criteria of Eq. (6) with a standard value of SD=0.
Assuming r 1 is the new x(t), imf 1 , imf 2 , …, and imf n are calculated successively until the last time series r n remains undecomposed. Thus, the original time series can be presented in terms of imf and a residual component r, as shown in Eq. (8).
where r n is the residual component representing the trend or mean of the original time series x(t).
The imf component is an oscillating function with different amplitudes and frequencies. Each imf component presents two characteristics [21]. In the data domain of each imf component, the numbers of the maximum values must be equal to the number of the zero crossings or present a difference of one at most. The average value of the envelope curves must always be equal to zero.
Taking the friction noise signal in test 4 as an example, the original signal is decomposed into 11 different imf components and a single r. Figure 4 shows the time-and frequency-domains of the original signal, imf 1 , imf 2 , imf 5 , r, and the reconstructed signal. Several peaks in the power spectra of the original signal and imf 1 and imf 2 components may be observed. By contrast, the power spectrum of imf 5 is relatively smooth. Therefore, the new signal is reconstructed by components from imf 5 to imf 11 and r. The power spectrum of the reconstructed signal is smoother than that of the original signal.    increases in the running-in process, remains relatively stable in the steady-state process, and then increases sharply in the increasing friction process. At this point, the specimens were damaged and the tests were ended.

Phase-space reconstruction
The phase space is a geometric space that reveals the states of a system. In general, a nonlinear dynamic system presents a very high phase-space dimension. However, in practice, data from the tests are considered to present a single-variable time series obtained from the interaction of different parameters in the system. Thus, the test data should be reconstructed into a high-dimensional space to gain more dynamic information. Takens [22] presented a method wherein a 1D chaotic time series is extended into a 3D or higher-dimensional phase space by phase space reconstruction. The aim of this method is to expose more information on the system hidden in the time series. Time difference method is usually utilized to reconstruct the phase space for 1D time series of x 1 , x 2 , x 3 , …, x n . A number selected from the original time series is as one of the components of the vector every τ times to construct a group of vectors, i.e, 2 ( [ , , , , ] i.e, where m is the embedding dimension (the dimension of the reconstructed phase space), X i is the vector of the reconstructed phase space, τ is the delay time, n is the length of the original time series, and N is the number of vectors in the reconstructed phase space.
Selecting an appropriate τ and m is very important for reconstructing the phase space. When τ is too small, x(t) and x(t+τ) cannot be independent of each other because their values are close to each other. When τ is too large, the relationship between x(t) and x(t+τ) becomes random, and the chaotic attractor cannot be accurately determined. In this regard, the autocorrelation function method is used to determine an optimal τ [23]. For a set of single-variable time series {x(i)}, the definition of the autocorrelation function method is The function constructed from τ and C(τ) is given by Eq. (10), where τ is the optimal delay time when C(τ) falls to the value of (1-1/e)·C(0). Taking the friction noise of test 4 as an example, Fig. 6 illustrates the relationship between τ and C(τ). Here, C(τ) falls with increasing τ. For example, the coordinates of the first point are (22, 0.6184), the τ of the time series is 22 s.
The precondition of choosing the m is m≥2d + 1 (d is the dimension of the dynamic system). For an infinite time series without noise, m is just larger than the smallest integer value of D. In an finite time series  with noise, however, m should be much larger than D. If m is too small, the attractor could undergo self-intersection because of folding. Thus, selecting a relatively large m is necessary in theory. Unfortunately, increases in m also increase the calculation burden for geometric invariants (e.g., D and the Lyapunov exponent) in practical applications. Moreover, the influences of noise and rounding errors significantly increase.
Therefore, an optimal m is selected by the saturated correlation dimension method [13]. The advantage of this method is that both D and m can be calculated at the same time. The dimension is calculated by means of the G-P algorithm, which was proposed by Grassberger and Procaccia [24]. The formula of the G-P algorithm is 1 1 , where H(·) is the Heaviside step function. That is when r approaches infinity, D is defined as 0 ln ( ) lim ln Reconstructing the phase space of a 1D time series is necessary prior to this calculation. The τ can be obtained by the autocorrelation function method described above, and the lnC(r)−lnr curves are plotted in double-logarithmic coordinates for each m. Ideallinear intervals are then selected from the curves and fitted by the least-squares method [25], the slopes of which reflect D. Then, m is plotted as the horizontal coordinate, and D is plotted as the vertical coordinate.
The slope values are the D corresponding to the m plotted in the curved figure. Finally, the optimal m and D are obtained when D becomes stable.
Taking the friction noise in test 4 as an example, D and m are calculated by the saturated correlation dimension method. Figure 7(a) shows that doublelogarithmic curves are obtained when the integers of m are from 11 to 30. Figure 7(b) displays the relationship between D and m obtained by fitting through the least-squares method. When m is less than 25, D increases with increasing m. When m is greater than or equal to 25, D remains stable over a small range. The optimal m and D of the friction noise in test 4 are 25 and 0.8977, respectively. Table 2 shows τ and the optimal m of friction signals from five tests.

Phase-space trajectory and attractor evolution
The τ and optimal m of friction noise signals were calculated by Takens' theorem [22] in Section 4.1. However, observing the evolution of phase-space trajectory visually in 3D space is impossible because of the high dimension of the evolution process. Therefore, the high-dimensional space is projected onto a 3D space by principal component analysis (PCA; Appendix A) [17,26]. Three main vectors selected from the reconstructed high-dimensional space are drawn in the 3D graph to observe the evolution of the phase-space trajectory expediently. The reconstructed vector X is given in Eq. (9). The inner product matrix Y of the vector matrix X is presented as follows: The three main eigenvalues λ 1 , λ 2 , and λ 3 are chosen from all of the eigenvalues of matrix Y and then calculated and arranged in descending order as λ 1 , λ 2, λ 3 , …, λ n (λ 1 ≥ λ 2 ≥ λ 3 ≥ … ≥ λ n ≥ 0). Then, a new N × 3 order vector matrix α is obtained by projecting the reconstructed m-dimensional vector matrix X to the three matrix directions V 1 , V 2 and V 3 , which are calculated according to λ 1 , λ 2 and λ 3 , respectively. 1 2 3 , , , , Each row of matrix α represents the coordinate of a point. A total of N points, which are drawn in 3D coordinates to present the 3D phase-space trajectory, are reconstructed.
Taking the friction noise signals in test 4 as an example, the τ and m of the friction noise are 22 and 25, respectively. The friction noise signal is divided into sections every 20,000 data points, and the points are drawn in 3D space by phase-space reconstruction. All of the phase trajectories are divided by continuous and non-overlapping signals. Because of space limitations, only some figures of the trajectories are selected in this work to present the evolution law of the friction noise signal. The phase trajectories of the friction noise are given in Fig. 8.
During the running-in process of the friction noise signal, the trajectory of the friction noise begins to converge, and the volume of the phase trajectory is very large. The radius of the trajectory is very large in the initial friction stage of 0-15 min, as shown in Fig. 8(a). In the 60-75 min stage ( Fig. 8(b)), the curvature radius of the trajectory gradually converges to a central point as the friction process continues. The attractor of the friction noise forms in this stage. In the 75-315 min stage, the trajectory of the friction noise converges to a smaller point, and the trajectory circles reciprocally. The curvature radius remains steady within a small range, as shown in Figs. 8(c)-8(g), and the attractor of the friction noise is stable. In the final process, the phase trajectory of the friction noise begins to diverge, and the curvature radius increases. Thus, the phase trajectory escapes from the space presented in Fig. 8(h), and the attractor of the friction noise disappears. The evolution of the phase-space trajectory of the friction noise can be defined as "convergence-stability-divergence", which corresponds to the "forming-keeping-disappearing" pattern of the evolution of the friction noise attractor. This pattern also corresponds to the complete friction process of the system. The friction noise first increases, and then declines to a stable value in the 0-75 min stage; this stage is considered the running-in process. The friction noise remains steady in the 75-315 min stage, which is also known as the steady-state process. Finally, the friction noise increases rapidly in the final process. During the complete friction process, the friction noise attractor gradually forms, then remains stable for a long period of time, and then finally disappears. The evolution of the attractor can be determined from the evolution of the phase trajectories. In the running-in process, the phase trajectory becomes convergent, and the attractor begins to form. During the stable-state process, the trajectory is maintained in a specific space, corresponding to the stable stage of the attractor. Finally, the phase trajectory diverges and escapes from the specific space, which corresponds to the disappearance of the attractor. Thus, the phase trajectory evolution of the friction noise is highly consistent with the evolution of the attractor. The attractor may be considered a running-in attractor because it forms in the running-in process.

Evolution of the correlation dimension
D is a type of fractal dimension and is sensitive to the time behavior of the system. To reflect the characteristics of the chaotic dynamic system, D is usually used to quantitatively describe the complexities of the chaotic signals on a small scale [17]. In the present case, D can be utilized to characterize the complexity of friction noise signals [17,27].
Coinciding with the time stages of the phase trajectory described in Section 4.2, the time series of the friction noise signal is divided into continuous and non-overlapping windows every 15 min (including 20,000 data points). Figure 9 illustrates the evolution of D of the friction noise in five tests as calculated by the saturated correlation dimension method discussed in Section 4.1.
In the attractor-forming process, D increases from a low value. In the attractor-keeping process, D remains stable. In the attractor-disappearing process, D declines. The evolution curve of D in the friction process corresponds to an "inverted bathtub curve", and can be generalized as "increasing-steadyingdeclining" consistent with the "forming-keepingdisappearing" formation process of the friction noise attractor and the "running-in, steady process, increasing friction" pattern of the complete wear process.

Evolution of the Lyapunov exponent
The basic characteristic of chaotic motion is that the system is extremely sensitive to the initial value. The trajectories generated by two initial values that are close to each other separate exponentially with increasing time. This phenomenon can be quantitatively described by the Lyapunov exponent. The Lyapunov exponent represents the average exponential rate of convergence or divergence between adjacent orbits in a phase space. When the Lyapunov exponent of a system is less than zero, the phase volume of the system is contractive in this direction. When the Lyapunov exponent is greater than zero, the phase volume is expansive and folding in this direction. The long-term behavior of an uncertain system is unpredictable. Therefore, the system is chaotic [28].
If a system has a chaotic attractor, it presents three features: (1) There is at least one positive Lyapunov exponent, (2) at least one of the exponents is zero, and (3) the sum of the exponential spectrum is negative.
In practice, the computational burden for all Lyapunov exponents is very large when the dimensions of the system are very high. Thus, determining chaotic characteristics through the largest Lyapunov exponent is appropriate [28].
The exact Lyapunov exponent of a general time series cannot be obtained according to the dynamic equation. The wolf reconstruction method [28] is thus used to calculate the maximum Lyapunov exponent of such a time series. The calculation process of this measure is as follows: The τ and m are respectively obtained by the autocorrelation function and the saturated correlation dimension methods. According to Takens' theorem, the new time series Y(t i ) = (x(t i ), x(t i+τ ), … , x(t i+(m-1)τ )), (i = 1, 2,…, N) is obtained by reconstructing the phase space of the original time series.
Assuming that Y(t 0 ) is the initial point in the phase space, Y 0 (t 0 ) is the point nearest to Y(t 0 ), and L 0 is the distance between these two points, the time evolution of these two points is tracked from t 0 until the distance is larger than ε at t 1 . Then, a point Y 1 (t 1 ) near the point Y(t 1 ) is determined to ensure that the distance between the two points satisfies the function.
To ensure that the angle between L 1 and L′ 1 is as small as possible, the process described above is repeated until the end time of Y(t). Assuming M is the iteration number, the largest Lyapunov exponent [29] is given as Eq. (19).
Coinciding with the time stages of the phase trajectory in Section 4.2, the time series of the friction noise signal is divided into continuous and nonoverlapping windows every 15 min (including 20,000 data points). The largest Lyapunov exponent of every stage is calculated by Eq. (19). Figure 10 illustrates changes in λ 1 of the friction noise over five tests. Because the λ 1 of each test is greater than zero, the friction noise signals are chaotic and the attractors of the friction noise signals can be considered chaotic attractors.

Conclusions
Experiments with different rotating speeds and loads were performed on a ring-on-disk tribometer, and friction noise signals produced during the friction process were extracted by a microphone PCB. The phase trajectories were projected onto a 3D space by PCA, and D and the largest Lyapunov exponents were calculated based on phase-space reconstruction. The following conclusions were confirmed.
(1) The chaotic attractor is a set of infinite points in phase space. The system state of chaotic motion always converges to a certain attractor in that phase space. Thus, the evolution of the phase trajectory can reflect the evolution of the chaotic attractor. The processes of convergence, stabilization and divergence characterize the evolution of the phase trajectories of friction noise signals, and the evolution of these phase trajectories corresponds to the evolution of the attractors via a pattern called "forming-keeping-disappearing". The formation process of the attractor can characterize the chaotic behavior of the friction system.
(2) D increases during the running-in process, remains relatively steady during the stable process, and then decreases during the friction increasing process. The change process of D conforms to the pattern of an "inverted bathtub curve". Changes in D are consistent with the evolution processes of chaotic attractors and the phase trajectories.
(3) The attractors of friction noise signals are also called running-in attractors because they form during the running-in process. In addition, the attractors are chaotic at the points where the positive largest Lyapunov exponents are obtained during the friction process.
(4) Friction noise includes information that can be utilized to characterize the dynamic behavior of a friction system. Therefore, the largest Lyaponov exponent can describe the chaotic characteristics of the friction process, and the phase trajectory and D of the friction noise can describe the friction process and evolution of attractors. The results of this study help reveal changes in wear states. As the evolutionary consistency of the chaotic characteristics of friction force and friction noise was not examined in this work, future research on the friction process will include this topic.
PCA is a statistical process that utilizes an orthogonal transformation to transform a set of observations of possibly correlated variables into a set of values of linearly uncorrelated variables. It is one of the simplest methods based on true eigenvector multivariate analyses and can reveal the inner structure of data. For a multivariate data set in a high-dimensional space, PCA can provide a lower-dimensional picture or a projection of this object when observed from its most informative viewpoint by applying only the first few principal components so that the dimension of the transformed data is decreased.
Mathematically, the definition of this transformation is a set of p-dimensional vectors of w (k) = (w 1 , ... , w p ) (k) that map each row vector x (i) of X to a new vector of principal component scores t (i) = (t 1 , ... , t m ) (i) , shown by where Σ represents the singular values of X, U represents the left singular vectors of X, and W represents the right singular vectors of X.
In terms of this factorization, the matrix X T X is Truncation of a matrix M or T by SVD produces a truncated matrix that is the nearest possible matrix of rank L.
Therefore, PCA can concentrate most of the signals into the first few principal components, which can be obtained by dimension reduction; later principal components may be affected by noise and disposed of without great loss of information.
Open Access: The articles published in this journal are distributed under the terms of the Creative Commons Attribution 4.0 International License | https://mc03.manuscriptcentral.com/friction (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.