On the boundedness of running-in attractors based on recurrence plot and recurrence qualification analysis

A feature parameter was proposed to quantitatively explore the boundedness of running-in attractors; its variation throughout the friction process was also investigated. The enclosing radius R was built with recurrence plots (RPs) and recurrence qualification analysis (RQA) by using the time delay embedding and phase space reconstruction. Additionally, the typology of RPs and the recurrence rate (RR) were investigated to verify the applicability of R in characterizing the friction process. Results showed that R is larger at the beginning, but exhibits a downward trend in the running-in friction process; R becomes smooth and trends to small steady values during the steady-state friction period, and finally shows an upward trend until failure occurs. The evolution of R, which corresponded with the typology of RPs and RR during friction process, can be used to quantitatively analyze the variation of the running-in attractors and friction state identifacation. Hence, R is a valid parameter, and the boundedness of running-in attractors can offer a new way for monitoring the friction state of tribological pairs.


Introduction
The friction coefficient (or "coefficient of friction", COF) generated from friction couples is an important information carrier of friction conditions, directly reflecting wear state and wear mode [1,2]. Owing to their strong dependence on a tribological system and their time varying nature, friction signals inevitably display nonlinear and complex characteristics [3]. The study of the COF based on nonlinear theory has gained much attention over the years, particularly with respect to the field of fractal and chaos [4].
The sequence signals extracted from friction systems, including friction force, friction coefficient, frictioninduced vibration, and friction temperature, have been proven to be chaotic, and the output phase spaces show obvious low dimension [5,6]. The phase trajectory, obtained by connecting the phase points in turn, was demonstrated to intuitively describe the time evolution process of a system [7,8]. The attractors are spontaneous ordered spatiotemporal structures of phase trajectory during the running-in process; thus, the attractors in tribo-system are also referred to as running-in attractors. As typical chaotic attractors, the running-in attractors have the properties of boundedness, fractal-dimensionality, and intrinsic randomness [9], which can be characterized by fractal and chaotic parameters, including fractal and multifractal parameters [10,11], Lyapunov exponent [12], entropy [13], and predictability [14]. Each measurement presents and emphasizes the attractors from a different perspective. It has been widely confirmed that both the fractal dimension and multifractal spectra play a vital role in characterizing the fractal structures, while the Lyapunov exponent contributes to measuring the speed of divergence or convergence of nearby orbits in phase space, and the entropy measures the complexity of the attractors. Additionally, Zhou et al. [15] had successfully introduced predictability as a quantitative statistical judgment of the intrinsic randomness degree in running-in attractors.
The study of running-in attractors focuses not only on their properties and characterization parameters, but also on their extensive use in wear condition identification. Zhu et al. [16] discovered that running-in attractors are shaped during the running-in process. Sun et al. [17] stated that the evolvement of the friction vibration attractors corresponds to changes in wear state from the running-in stage to the stable wear stage. Zhou et al. [18] indicated that the evolution of friction temperature attractors abided by the law of "forming-keeping-disappearing", which was consistent with the stages of "running-in, steadystate, and increasing friction". Liu et al. [19] analyzed the cross correlation between tangential and normal friction vibration, and suggested that the chaotic attractors of both signals can be used to describe variations in the running-in process. Jiang et al. [20] developed a reliable method for fault diagnosis for gear boxes based on the chaotic degree of frictioninduced vibration attractors; the results showed that the chaotic attractor was maintained in a steady state, and the chaotic degree was related to fault severity.
The above-mentioned information and references are only a fraction of existing research papers that discuss running-in attractors, and there remains some unresolved issues. Firstly, as noted in Ref. [16, 18−20], the friction states can be recognized by the evolution of the phase trajectory of running-in attractors. However, there is no quantitative description of the degree of convergence of the running-in attractors, thus the existing studies on wear condition identification by phase space trajectory merely rest on qualitative analysis [21,22]. Additionally, effective characterization parameters to appraise the boundedness of the running-in attractors are lacking. Zhang et al. [23] indicated that the bounds of a domain containing all compact invariant sets obtained in many cases can be used not only for theoretical studies of chaotic attractors, but also for estimating for the Haudorff dimension and/or for further numerical researches. Furthermore, it has been ascertained that there are obvious shortcomings in the characterization of rather short-time nonstationary sequences with the method of entropy, fractal dimension, and the maximum Lyapunov exponent [24,25].
The primary goal for this paper is to put forward a feature parameter to quantitatively characterize the boundedness of running-in attractors, and to further investigate the evolvement of the nonlinear dynamic behavior of tribological systems by describing the variation in the proposed parameter. Moreover, recurrence plots (RPs) and recurrence quantification analysis (RQA) are explored to verify the proposed parameter. This paper is organized as follows. In Section 2, the friction experiments are described, and the COF signals are acquired throughout the friction process. In Section 3, measurement of the enclosing radius R based on RPs and RQA is introduced to define the boundedness of running-in attractors. In Section 4, the morphological evolution of RPs, and the evolution of feature parameters are further studied. Finally, the major conclusions of the paper are summarized in Section 5.

Apparatus
A self-regulating rotational motion tribometer illustrated in Fig. 1(a) was applied to conduct friction experiments. It can be clearly seen that the upper ring sample ( Fig. 1(b)) was fixed into the ring holder by a locking bolt and locating pin, while a center hole and a mill groove were machined on the disc holder to locate the lower disc sample ( Fig. 1(c)). Moreover, the wedge-ended pin was used to ensure relative stability of the disc sample. Rotary motion of the main shaft, driven by a motor, was carried through the ring holder to the upper sample, causing relative sliding between the contact surfaces. A frequency converter was applied to realize rotational speed adjustment, and the normal load was adjusted by weight. The COF was then calculated via the collected friction torque signal, which was measured using the static torque sensor mounted on the disc holder bottom, and then uploaded to the computer via the signal collecting device. With μ as the COF and M as friction torque, the relationship is given via the following formula [26]: where R 1 and R 2 are the outer and inner radius of the ring specimen, respectively. P is the normal load applied on the contact surface.

Test samples and test conditions
The lower disc sample was made up of AISI 1045 steel with an HB194 hardness and was cut into Φ 46 × 12 mm segments. The upper ring sample was made of AISI 52100 steel that was hardened up to 698 HB. The ring sample had an outer diameter of 34 mm and an inner diameter of 24 mm. Thus, the nominal contact area was 455.53 mm 2 , and the equivalent radius was 14.5 mm.
The initial surface roughness of the ring and disc were 0.036-0.040 μm and 5.60-5.90 μm, respectively. It is noted that both rings and discs were carefully cleaned with anhydrous ethanol via ultrasonic bath prior to being used in the friction experiments. Four tests were carried out under different normal loads and velocities, and the experimental conditions were summarized in Table 1. The samples were lubricated with CD 15W-40 type machine oil with a viscosity of 12.5-16.3 at 100 °C , with the flash point temperature higher than 215 °C and freezing point temperature lower than -20 °C . Prior to tests, a volume of 0.2 mL lubricating oil was added to the contact surfaces and no more oil was supplied during

Test results
Any processed signal inevitably contains an abundance of noise, which can not only degrade the accuracy of parameters as-calculated, but also obscure or even destroy the inherent principles of the signal itself. Empirical mode decomposition (EMD) is widely known to be capable of eliminating the noise effects, and was introduced in the study for the de-noising. The basic functions of the EMD method are derived from the signal itself, so the method can be used to analyze nonstationary and nonlinear processes [27]. By using EMD, noisy signals can be decomposed into a number of intrinsic mode functions (imfs) and a residue r. Each imf has only one frequency, and can reflect the intrinsic and real information of the signal [28]. Taking Test 3 as an example, the temporal waveform of the original COF signals and its power spectral density (PSD) spectra are respectively shown in the top left-hand-side and right-hand-side diagrams of Fig. 2. It can be clearly observed that the raw COF signals contained two frequency components before decomposition. Therefore, the first two decomposed imfs (imf1 and imf2) and their PSD spectra are shown in the successive left-hand-side and right-hand-side diagrams, respectively. Obviously, the waveforms of imf1 and imf2, and their PSD spectra displayed in Figs. 2(c)-2(f) show that they are nearly monocomponent, whilst the frequency of the first two imfs coincides with the dominant frequencies of the original signals, which is approximately 110.2 Hz and 79.2 Hz, separately. Thus, it can be confirmed that these two imfs mainly correspond to the noise components.
The de-noised COF signals throughout friction process of four tests are shown in Fig. 3. Although the wear lives in four tests are not the same on account of different working conditions, the COF signals evolve in a similar way, which can be mainly divided in three stages: the running-in, steady-state, and increasing friction stages. The COF signals first increased, quickly reaching a peak value in the running-in friction process ( Fig. 3(b)(i)), then gradually decreased and remained in a stable value in steady-state friction process ( Fig. 3(b)(ii)), and the COF sharply increased in the increasing process ( Fig. 3(b)(iii)). Additionally, there were fluctuations in signals at mean value, especially in the running-in and increasing friction stages.

Phase space reconstruction
Reconstructing the friction signal time series into an appropriate high dimension phase space is an effective means to investigate the chaotic characteristics of the friction process in a tribological system. According to the time delay embedding technique proposed by Takens, the method of delays can be used to embed a scalar time series {x i } (i=1,2,3,…,N) into an m-dimension space as follows [29]: The phase space can be expressed as follows: where N is the number of the phase space vector and N=n-(m-1)τ, τ is the delay time and m is the embedding dimension.
Takens' theorem assumes that the time delay τ can be chosen almost arbitrarily in the case of an infinite noise-free data set [30]. However, the real data sets are finite and noisy; thus, the selection of embedding parameters plays a significant role in the construction of the phase space. In the presented work, τ and m are established by applying the mutual information (MI) method and the false nearest neighbors (FNN) method, separately. As suggested in Refs. [31,32], τ can be determined as the first minimum of the MI function, and according to Refs. [33,34], the phase space is entirely opened when the percentage of the FNN point decreases to 5%, and the optimal embedding dimension m is computed.

Recurrence plot and recurrence rate RR
A recurrence plot (RP), a visualization tool introduced by Eckmann et al. [35], is presented to visualize recurrences in multi-dimensional phase space. Concretely, it is a two dimensional graphical representation of the trajectory of the dynamical system in the form of a binary recurrence matrix R, and the entry R i,j is expressed as follows: where R i,j is an element of the recurrence matrix; N is the number of the vector {X i } (i=1,2,3,…,N); ε is the threshold distance; the notation ║·║ is calculated as a Euclidean norm in the phase space, and Θ(·) is the Heaviside function, defined as follows: Owing to unnecessary prior assumptions on statistical properties, the RPs, along with RQA, are applicable to the characterization of nonlinear time series and have gained an increasing amount of attention in a variety of fields [36,37]. The recurrence for a time series is when a point on the trajectory repeats itself, which means that a point is close enough to another one within a suitably selected interval of error [38]. The basic measure of the RQA is the recurrence rate (RR) [39], which measures the percentage of darkened pixels in an RP: RR is a measure of the density of recurrence points and corresponds to the correlation integral. However, the calculation of correlation dimension needs a large amount of data points, and the determination of the scaling region is relatively subjective. We can avoid the constraints with RR. By employing the definition of the R i,j in Eq. (4), we can directly relate RR to ε: An RR value, ranging from 0 to 1, is a measure based on overall probability that a certain state recurs in an RP. Given a specified ε, a large RR indicates an intensive distribution of phase points and a high degree of convergence. Furthermore, a trajectory with a larger RR has a stronger space-filling capacity.

Feature parameter R of chaotic attractors
In order to extract a quantitative characterization of the size range where running-in attractors exist, the recurrence matrix and RPs are obtained based on phase space reconstruction, and the relationship between RR and ε is defined as in Eq. (7) is calculated. The cutoff distance ε defines a hypersphere centered at phase points, every point in the phase space is a neighbor of every other point when ε is large enough, and the phase space will be included in a hypersphere of radius ε, which can be defined as the enclosing radius R. R can be computed by the plot of RR and ε, and the eigenvalue corresponds to the first abscissa when RR = 1. The phase points can be completely enclosed when ε is greater than R.
We take the typical Lorenz chaotic attractor, which is defined by Eq. (8), as an example to illustrate the calculation of R. The differential equations are solved numerically by the four-order Runge-Kutta method with a step of 0.01 and iterating time is 150 s, after discarding the initial transient points to allow the trajectory to fall to the attractor. A file containing 4000 data points is obtained as a time series.
The determination procedure of embedding parameters of the Lorenze time-series ( Fig. 4(a)) are shown in Figs. 4 (b)-4(c), τ corresponds to the first minimum of the MI function, and the percentage of FNN stabilizes after m = 4, as indicated by the filled dot. Therefore, for the given time series, τ= 15 and m = 4 are used for the phase space reconstruction and the following analysis. The RPs of the Lorenz time series dependent on the cutoff distance ε are shown in Fig. 5.
As shown in Fig. 5, a fixed ε results in a symmetric RP, and the RPs exhibit characteristic small-scale patterns that are caused by different ε values. When ε is relatively small, there are few recurrence points, Fig. 4 One-dimensional time series y of the Lorenz system and embedding parameters. and the RP is sparse, with an increase in ε, the pairs of recurrence point increase, and the black area is larger than the white area within an RP. Figure 6 illustrates how the ε can determine RR. RR was calculated by a fixed ε and using a Euclidean distance. From results shown in Fig. 6, RR first increases and then tends to stabilize with increasing ε, and R of the given theoretical time series is 61.72. It is important to note that the magnitude of the parameter is not related to whether the attractor is chaotic, but is a quantitative representation of the volume and degree of convergence of the attractors. A lower magnitude of R indicates that in a smaller volume where the attractors exist and a higher degree of convergence of phase points.

Recurrence plots of COF signals
As noted above, such recurrence of a state at time i at a different time j is pictured within an RP with black and white dots, where black dots mark a recurrence. It has been generally confirmed that the visual inspection of RPs reveals four typical large-scale patterns (the typology): homogeneous, periodic, drift, and disrupted [40]. The typology of the RPs can qualitatively give hints about specific dynamic behaviors of a system. The first two patterns correspond to stationary systems with short relaxation times and periodic (or quasi-periodic) systems, respectively. On the contrary, non-stationary systems with slowly varying parameters have RPs paled away from the LOI (line of identity with an angle of π/4) and can be classified as drift, while abrupt changes or extreme events in dynamics cause white areas and lead to a disrupted RP.
The embedding parameters of the four tests were established by applying the method introduced in Section 3, and the results are summarized in Table 2. The recurrence matrices of the COF signals were calculated by applying a threshold such that RR was 50% and using the Euclidean norm, the RPs are then obtained via the MATLAB CRP Toolbox [ 41 ] . For shorter computational time and smaller memory requirements, the COF signals were divided into non-overlapping subsequences with an equal length of 23785 data points, and only the initial 3785 points were calculated at every subsequence. Owing to space limitations, this paper enumerates only the 16 profiles in Test 2, as displayed in Fig. 7.  An overview of the RPs' evolvement in the running-in attractors throughout the friction process is illustrated in Fig. 7, and the macroscopic patterns of RPs depict an evolution of "disrupted-drifteddisrupted". It can be observed in Figs. 7(a)-7(d) that large black and white areas occur at the corner of the RPs, indicating the friction system is unstable during the period. The disrupted typology corresponds to a violent fluctuation and a sparse distribution of phase points. In Figs. 7(e)-7(m), there is a fine structure of vertical (horizontal) lines, essentially demonstrating that the phase space trajectory does not change or changes very slowly. It can be recognized as typical behavior of laminar states [42]. Interestingly, there is a checkerboard structure in Fig. 7(n); it seems plausible that on this time-scale the underlying states strongly deviate from the previous "laminar" states. After wearing for 238 min, the COF signal changes immensely and a large fluctuation starts to appear. The RPs shown in Figs. 7(o)-7(p) are characterized with large black regions concentrated at the lower left corner and can be regarded as disrupted patterns.
The RPs evolve from a disrupted pattern to a drifted pattern during the running-in process, maintain the drifted pattern in the steady-state process, and finally return to a disrupted pattern in the increasing friction process. Although the amplitude analysis of the friction coefficient was more likely not accurately revealed, the qualitative features can be extracted for comparison.

Feature parameters of running-in attractors
In order to illustrate the applicability of the proposed running-in attractors quantitative measurements of COF, R and RR are computed for each signal, and the evolution of R and RR are displayed in Fig. 8.
As shown in Fig. 8, in different friction stages, the characteristic parameters R and RR show temporal variation and preferable regularity. Within the running-in friction process, R shows a downward trend; during the steady-state friction period, R becomes smooth and trends towards the small steady values; and R shows an upward trend until failure occurs. On the contrary, RR increases at first, then maintains stability, and finally increases. It is shown that both the R and RR values computed from different tests evolve in a similar way, and the evolution is synchronous if both eigenvalues are extracted from one friction coefficient signal. The results of R and   Fig. 7.
The Pearson correlation r is adapted as the statistic for measuring the association of the changes in two parameters. As shown in Table 3, the Pearson correlation r of R and RR obtained by the four tests are negative, and the absolute values are all greater than 0.76, indicating there is a fine negative correlation between the two parameters.

Discussion
The phase trajectory is an effective tool for describing the time evolution process of a system intuitively; thus, it is perfectly suitable for analyzing tribological issues. According to the evolution of the phase trajectory of friction signals, the friction process can be interpreted as "self-organization, chaos, and system-instability", which correspond to the patterns of "running-in, steady-state, and increasing" across the friction process.
Owing to the rather rough initial machined surfaces, the contact pressure was large and the friction was severe at the very beginning of the experiments. The COF significantly fluctuated, indicating a sparse distribution of phase points and a high degree of divergence. The typology of RPs could be interpreted as disrupted, and the bounds of friction coefficient attractors were rather large. R was rather large, while RR was small.
Within the running-in friction period, asperities grew increasingly abundant and came into contact; the friction coefficient grew stable with decreasing real contact pressure. The distribution of phase points becomes relatively concentrated, and the bound of the phase trajectory gradually converged. Hence, the running-in process can be regarded as the formation process of running-in attractors. This results in a uniformly distributed RP, a decrease in the shrinking of R, and an increase in RR.
In the steady-state friction stage, the friction system was well lubricated, and chiefly characterized by steady surface contact with low and stable heat generation. When wear failure occurred, the COF increased and fluctuated sharply as a consequence of the generation of wear particles, as well as the deterioration of lubricant and worn surfaces. The phase trajectories gradually diverged and finally escaped the limited space in the increasing friction process, causing a disrupted RP. Moreover, R increased and RR decreased dramatically.
Furthermore, for the same pressure, the amplitude and the fluctuation of R both decreased with increasing velocity. Meanwhile, for the same velocity, the amplitude of R decreased while the fluctuation increased with an increase in normal pressure. In the process of sliding friction, the actual contact and interaction of surface asperities can be dramatically reduced by the increasing velocity, leading to a stable friction coefficient. With increasing pressure, the actual contact surface and the number of interacted asperities increase, causing a smaller friction coefficient amplitude. At the same time, it can be stated that a higher pressure will result in an increasing mechanical engagement between the asperities, as well as a fluctuation in the friction coefficient.

Conclusions
The friction experiments were performed by sliding a ring (AISI 52100) against a stationary disc (AISI 1045) under lubricated conditions. The evolution in running-in attractors was investigated by defining a feature parameter R extracted by RP and RQA. The dynamic evolution of friction systems was intuitively analyzed by RPs and characteristic parameters R and RR. The major conclusions can be summarized as follows: (1) The enclosing radius R, defined by RP and RQA, is adapted to describe the boundedness of running-in attractors. It is an effective way to reflect the behavior of COF during friction process. A small value of R indicates an intensive distribution of phase points and a high degree of convergence, while the convergence of trajectory indicates a decrease of energy dissipation.
(2) R goes through a rapid downward path during the running-in period, a stable fluctuation in the steady-state friction stage, and an uptrend in the increasing friction stage, following the same trend as the COF in time domain. The identification of a friction state with R is consistent with that using the COF signals; therefore, R can be used to monitor the friction state of tribopairs.
(3) The results of RPs and RQA indicate that the non-linear dynamic behavior of the friction process can be described by the pattern of RPs and RQA measures. The drifted typology of RPs, the decrease in R, and the increase in RR occur when the running-in process is terminated. Additionally, the disrupted typology of RPs, the increase in R, and the decrease in RR indicate the failure occurs.
(4) The typology of RPs and chaotic parameters of COF signals during different friction tests evolve in a similar way, which can be used for wear condition identification. Furthermore, they synchronously evolve when the signal is extracted from a friction system.