Nonlinear dynamic characteristics of load time series in rock cutting

2.1. Phase space reconstruction

The irregular contact between the cutting mechanism and rock usually contributes to the complexity, nonlinearity, ambiguity, and dissipativeness of a cutting system and to the nondeterminacy of a cutting load time series. However, the complex dynamic characteristics of a cutting system are difficult to verify in a low-dimensional coordinate systems. To investigate such dynamic systems, the time series must undergo phase reconstruction. However, the cutting load signal can only be collected usually during rock cutting, and thus, obtaining reconstructed phase space along with all the derivatives of the cutting load signal can yield plenty of errors. A time series is determined by the self-relevant components defined in the Takens theorem [8], and the information on its components is hidden within the time series itself. Therefore, one or several equivalent phase spaces of the time series can be built using the time delay method to recover the information on the cutting system. Several methods have been developed recently for phase-space reconstruction, such as uniform and non-uniform embedding procedures [9-10]. The delay-coordinate method for phase-space reconstruction of 1D time series ($x\left(i\right)$, $i=\mathrm{1,2},\dots ,n$) has been widely used in different fields by Packard [11], and the $m$ – dimensional reconstructed vectors:

where $m$ represents the embedded dimension of restructured phase space; $\tau$ stands for the delay time; $N$ is the points number in the phase space, $N=n-(m-1)\tau$.

2.2. Denoising based on wavelet analysis

As a result of the influence of measurement devices and the accuracy of signal acquisition, the desired signal of the cutting load is often disturbed by noise. Consequently, the internal dynamics of the cutting system may be masked by the unpredictability and destructive effect of noise. Therefore, noise must be reduced to reveal the internal characteristics of the cutting load. Numerous methods have been developed recently for signal denoising, such as filtering, shadow theorem, singular spectrum analysis, and wavelet analysis [12-13]. Among these methods, wavelet analysis is considered as the most suitable for signal denoising [14]. Denoising analysis is generally conducted in either the time or frequency domain, and thus, wavelet analysis can be conducted in both domains, which helps in processing random and non-stationary time series to achieve a high signal-to-noise ratio. For a discrete signal $x\left(i\right)$ with a family of shifted and scaled wavelets associated with $\phi \left(t\right)$, the discrete wavelet transform is given as [15]:

The mother wavelet in Equation (2) is discretely dilated by the scale parameter $a_o^g$ and is translated by the translation parameter $i{b}_o{a}_o^g$, where $a_o$ and $b_o$ are fixed values, with $a_o>1$ and $b_o>0$, whereas $g$ and $i$ are positive integer values [16].

2.3. Delay time and embedded dimension

Delay time is an important parameter for phase-space reconstruction. The trajectories in reconstructed phase space gather at the same place if the delay time was too short because of the influence of the acquisition instrument and interference noise on the cutting load time series. Moreover, the shortest time interval for capturing the cutting load information cannot be covered. If the delay time is too long, the information cannot be captured because of the violent fluctuations in the cutting system. Several methods, such as the autocorrelation and mutual information methods can be used to estimate the delay time. The autocorrelation method is preferred for its simplicity and maturity in estimating delay time. For a 1D discrete time series, the autocorrelation function of delay time $j\tau$ can be expressed as [17]:

Given $j$, ratio of the autocorrelation function with delay time $j\tau$ to that of the initial value is solved. $\tau$ can be universally considered as delay time when $R_{xx}\left(j\tau \right)/R_{xx}\left(j\right)$ decreases to $1-1/e$ or if it reaches the first minimal value [18].

Embedded dimension is another important parameter for phase-space reconstruction. The trajectories in the reconstructed phase space fold if the embedded dimension is too small. If the embedded dimension is too large, then the restructured phase space increases the calculated amount, and noise may be dominant. Several methods can be used to estimate embedded dimension such as the correlation integral, nearest neighbour and odd value decomposition methods. This study uses the correlation integral method because of the high density of points in the phase space. Given $r>$0 as threshold, the number of points according to $‖x_i-x_j‖<r$ is calculated. The ratio of the calculated points to that of the total points is the correlation integral, which can be expressed as [19]:

where $H\left(*\right)$ represents the Heavside function $H\left(>0\right)=1$, $H\left(\le 0\right)=0$.

The threshold $r$ is set to an appropriate value, and the correlation dimension $D$ of the time series can be calculated approximately as:

2.4. The Lyapunov exponent

The power spectrum aims to reveal the frequency characteristics of the time series in the time domain or frequency domain as well as to distinguish periodic, quasi-periodic, and random signals. Despite its intuitionism and easiness, the power spectrum cannot readily distinguish the large periodic dynamic system or chaos system. The Lyapunov exponent is used to describe the diffusion or convergence rate of the phase-space trajectory of chaotic time series with an exponential law. The chaos system is confirmed by a positive Lyapunov exponent. This exponent is also widely used to identify the chaotic characteristic of time series. The Wolf, Jacobian, small data methods can be used to calculate the maximum Lyapunov exponent at present. The small data method has several advantages, including its high-speed calculation and low sensibility to embedded dimension, delay time and data amount. For phase-space reconstruction, the nearest neighbor of each point can be found in phase space by using the following equation [20, 21]:

where $t$ is equal to $\mathrm{1,2},\dots ,N$ and $t\ne \widehat{t}$; $p$ is the average period of the cutting load time series, that can be set to the reciprocal of the energy spectrum frequency; $X\left(t\right)$ is the state point of the reconstructed phase space; $d_t\left(0\right)$ is the distance between one point and the nearest point during the initial time.

Suppose that an exponent diverging rate $\lambda$ is present between reference the point $X\left(t\right)$ and its nearest point $X\left(\widehat{t}\right)$ in the basic orbit, then it can be calculated as follows:

where $\mathrm{\Delta }t$ is the sample period; $d_t\left(i\right)$ is the distance of the nearest points of $t$ group in the basic orbit with the elapsed time $i\mathrm{\Delta }t$; $〈\text{ln}{d}_t\left(i\right)〉$ is the average value of all about $t$.

2.5. Fractal dimension

The strange attractor generally has a self-similar structure in the phase space, and the chaotic time series also has a statistical similarity to the fractal time scales [22]. Fractal dimension is an important parameter in describing the complexity of a system. The load time series can be regarded as an open curve in two dimensions, and its outline has a fractal feature. The fractal dimension can be used to reflect the geometric feature to analyze the complexity and irregularity levels. Topological, capacity, similarity, and box dimensions are considered as fractal dimensions. The box dimension method is a simple and mature method technique for estimating the fractal dimension. The number $N\left(l\right)$ of boxes with side length $l$ covers the curve of the time series $F$, which follows the fractal distribution with a different side length, as follows [23]:

According to Equation (8), fractal dimension is equal to the slope of $\text{log}2\left(N\right(l\left)\right)$ to $\text{log}2\left(l\right)$.

3.1. Cutting load signal acquisition

The signal acquisition system of rock cutting is shown in Fig. 1. The system mainly included a cutting device, a fixed beam for attaching strain gages, a signal acquisition instrument, a computer to process signals, and other assistive devices. In the cutting device, the cutter would act on the rock, and its cutting speed and depth could be adjusted. The strain gage was attached to the fixed beam and not directly to the cutter to adjust cutting parameters conveniently. The INV306U instrument, which was developed by Beijing Dongfang Vibration And Noise Technology Research Institute, was used to collect the signals to achieve voltage signal acquisition and amplification by using different types of strain gage bridges. The sample frequency was set to 1024 Hz to obtain more comprehensive information. The signal was processed through a personal computer with intelligent data acquisition and software DASP-V10 software. A full-bridge circuit of strain gages was used to collect voltage signal [24]. The changes in voltage signal were eliminated through beam bending when the beam torque was used to calculate the cutting load.

Fig. 1The signal acquisition system: (a) a cutting device; (b) a fixed beam and strain gage bridge; (c) a signal acquisition instrument; (d) a computer to process signals

3.2. Experimental phenomena

In rock cutting process, the time series of a typical load signal was collected as shown in Fig. 2. The extrusion of the cutter on the rock at the initial stage produced a crushing zone, and the cutting load increased nonlinearly along with penetration displacement. The cutting load continued to increase after the crushing zone was formed, and several small peaks exhibited that a number of local breakings would occur before the peak cutting load could be achieved. The rock cutting load instantaneously decreased to a small value after the peak cutting force was achieved, thus indicating that a large rock fragment was formed because of the extrusion of the cutter on the rock. The aforementioned phenomena are consistent with the rock cutting theory or the experimental results [1, 25-26]. The formation of crushing zone, local breakings, and the avalanche of large fragments alternately occurred, which contributed to the appearance of several smaller pecks between two adjacent large peaks. The aperiodicity of the cutting load amplitude and wavelength confirmed the discreteness and irregularity of rock cutting and breaking.

Fig. 2Cutting load time series

4.1. Chaos characteristics

Power spectrum analysis can be used to understand the characteristics of a time series, such as periodic, quasi-periodic and non-periodic signals. The cutting load time series was translated by fast Fourier transform in the frequency domain [27], and its power spectrum is shown in Fig. 3. The power spectrum of the periodic signal exhibited discreteness; however, that of the cutting load exhibited a continuous broadband, non-dominant frequency, and a wide-band noise. Therefore, the cutting load time series exhibited chaotic characteristics. Given that such features were obtained outwardly, phase-space reconstruction and the maximum Lyapunov exponent were still necessary to identify the characteristics of cutting system.

The phase space of the time series was restructured using Equation (1), as shown in Fig. 4. The trajectories were disorganized and unsystematic in the phase space, and exhibited pseudo-randomness. Therefore, the interference noise should be reduced to prevent the time series from being influenced by the unpredictability and destructiveness of noise.

Fig. 3Power spectrum of cutting load

Fig. 4Reconstructed phase space

In the present study, the principle of denoising was adopted based on a multi-resolution signal decomposition analysis. The rock cutting signal was processed through multi-resolution signal decomposition, and was decomposed up to level 8 by using “discrete Meyer” as a mother wavelet. After inspecting the components, the signal was denoised by using detail component 7 and 8. The cutting load signal of denoising is shown in Fig. 5 whereas the reconstructed phase space is shown in Fig. 6. The trajectories of the reconstructed phase space were free from repetition, and enfoldment, and were concentrated at a certain region. Moreover, the reconstructed phase space of denoising exhibited an unambiguous strange attractor, which confirmed the chaotic dynamic characteristic of the cutting load. The denoising of the time series by wavelet analysis was shown to improve the deterministic components of the cutting system.

Fig. 5Cutting load signal of denoising

The ratio of the autocorrelation function with delay time $j\tau$ to that of the initial value could be obtained using Equation (3), as shown in Fig. 7. The delay time was equal to 10 when $R_{xx}\left(j\tau \right)/R_{xx}\left(j\right)$ reached the first minimal value, and to 100 when $R_{xx}\left(j\tau \right)/R_{xx}\left(j\right)$ decreased to $1-1/e$. In the presented study, 100 was selected as the delay time to avoid the loss of dynamic information of the cutting system.

Fig. 6Phase space of denoising signal

Fig. 7Rxx(jτ)/Rxx(j)

The correlation dimension saturated to a certain value with increasing embedded dimension if there was a strange attractor in the restricted phase space. The saturation value is often regarded as the correlation dimension of the time series. The variation of the correlation dimension with increasing embedded dimension is shown in Fig. 8. According to Takens theorem, the embedded dimension of the cutting system should be equal to 8.

The Lyapunov exponent of the time series can be estimated according to the slopes of the straight lines based on Equation (7). The maximum Lyapunov exponent was equal to the slope of a linear area of $y\left(i\right)~i$ according to least-square method. According to $y\left(i\right)~i$ and $y\left(i\right)-y(i-1)~i$ (Fig. 9), $y\left(i\right)-y(i-1)~i$ were equal basically equal when $i$ was within the range of 6 to 12. The slope of 0.065 within this range was regarded as the maximum Lyapunov exponent of the cutting load time series by the linear least-square regression method, which also demonstrated the chaotic dynamic characteristic of the cutting system.

Fig. 8Correlation dimension

Fig. 9Lyapunov exponent

4.2. Factors that influence fractal dimension

4.2.1. Cutting depth

According to Equation (8), the double logarithmic curves with the cutting depths of 15 mm and 20 mm are shown in Fig. 10. The slopes of these curves were 1.57 and 1.46 respectively, via the linear least-square regression method. The correlation coefficients of the fitting straight lines were all higher than 0.99, which reflected the fractal feature of the cutting load time series. Therefore, cutting depth had no influence on the fractal distribution of the load signal. A difference was observed in the fractal dimension when the cutting depth was different. The fractal dimension decreased with increasing cutting depth, thus indicating an increase in the self-similarity of cutting load when the cutting depth was reduced. Based on the size and number of rock fragments, a small fractal dimension produced large size and small number of rock fragments, thus indicating that rock was not seriously broken. By contrast, a big fractal dimension produced large number and small size of rock fragments. Therefore, the fractal dimension of the cutting load time series could be used to evaluate the extent of rock cutting or crushing.

Fig. 10Fractal dimension of cutting load with different cutting depth

(a) 15 mm

(b) 20 mm

4.2.2. Cutting velocity

The double logarithmic curves with the cutting velocity of 0.5 m/min and 4 m/min are shown in Fig. 11. The correlation coefficients were all higher than 0.99, thus indicating that the cutting velocity did not affect the fractal feature of the time series. The fractal dimension was equal to 1.55 at the cutting velocity of 0.5 m/min, and to 1.48 at the cutting velocity was 4 m/min. The time series with low cutting velocity had high self-similarity, and was capable of filling the phase space of the cutting system. Based on energy dissipation, a bigg fractal dimension with a low cutting velocity highly dissipated the energy of cutting system. With regard to size distribution, a low cutting velocity produced small fragments and big fracture areas, thus indicating that the energy required to generate a new fracture surface is big based on Griffith theory [28]. A low velocity could break down big fragments more than twice, and a high velocity might cause big fragments fly off from the motion path of cutter, thus preventing them from being broken down further. Therefore, the fractal dimension of the time series with a high velocity was smaller than that of the time series with a low velocity. Based on theory of rock crushing, the extent of rock crushing is high when the impact velocity is also high. However, the opposite is true regard to cutting velocity, as mentioned earlier. The present phenomenon could probably be attributed to rock fragments being crushed more than twice by the cutter during rock cutting. Therefore, the influence of cutting velocity on rock crushing and cutting performance should be investigated further.

Fig. 11Fractal dimension of cutting load with different cutting velocity

(a) 0.5 m/min

(b) 4 m/min

4.2.3. Types of assisted waterjet

High-pressure waterjet has been widely used in rock cutting, however, the fracture mechanism of the involved process remains unclear because of the opacity and damage instantaneity of the rock. To determine the rock fracture mechanism with different types of assisted waterjet (Fig. 1), this study investigates the fractal dimension of the cutting load time series. The double logarithmic curves under the waterjet front of the cutter or behind of the cutter are shown in Fig. 12. The correlation coefficients were higher than 0.99, thus indicating that such types of the assisted waterjet had an insignificant on fractal feature. The fractal dimension of the cutting load changed remarkably when the fracture mechanism of the rock was changed [3]. However, the difference in fractal dimension of the two different types was small, thus indicating that the fracture mechanism of rock cutting was unchanged by these types of waterjet. The waterjet in front of the cutter could decrease cutting load because of the the damage on the caused by the impact of the waterjet [29], however, the waterjet is used to clear fragments and to lubricate the cutter when it is behind of the cutter.

Fig. 12Fractal dimension of cutting load with different assisted waterjet types

(a) Waterjet front of the cutter

(b) Waterjet behind of the cutter