Recurrence analysis of regular and chaotic motions of a superelastic shape memory oscillator

The recurrence analysis is a promising tool for diagnostics of periodic and chaotic solutions, as well as identifying bifurcations. This paper deals with the application of this analysis for the first time to identify regular and non-regular motions of a superelastic shape memory alloy oscillator. The numerical analyses show that the method is capable of distinguishing periodic and chaotic trajectories. Recurrence quantities are applied, showing that different approaches are possible to establish the distinction between periodic and chaotic signals. Basically, recurrence entropy, trapping time, and characteristic recurrence time are considered.


Introduction
Shape Memory Alloys (SMAs) exhibit a complex thermomechanical behaviour determined by the occurrence of phase transformation between solid phases, namely austenite and martensite, at the microscopic scale [1,2]. This produces several functional properties, like superelasticity (or pseudoelasticity) and shape memory effects, which are largely utilized nowadays, for several applications in various fields. In particular, the superelastic behaviour is characterized by hysteretic loops without residual displacements ( Figure 1).
The force-displacement curves ( Figure 1a) are characterized by two plateaus associated with the occurrence of the Forward A→M transformation (FwT) and Reverse M→A transformation (RvT). The hysteresis loops can, in general, have different shape in tension and compression. Figure 1b also shows the temperature dependence of the four transformation forces that characterize the hysteresis loop, as well as the forcetemperature path that corresponds to the pseudoelastic loop in (a). In particular, FwT is associated with temperature rises, while RvT to temperature decreases. The intersections of the various lines with the temperature axis in (b), show four transformation temperatures -MS, Mf, AS, Af -that characterize a specific SMD. By and large, depending on the specific SMA, each of the main properties that characterize the pseudoelastic loops can assume different values in tension and compression.
When employed in dynamical systems, superelastic SMA elements may induce non-regular chaotic responses, as shown for example in [3][4][5][6][7][8][9][10][11]. The eventual presence of such chaotic motions may be an important aspect considered in the design of application and, for this reason, it is very important to be able to detect them. Due to the hysteretic behaviour induced by the non-smoothness of underlying differential equations, the most typical tool for chaos diagnosis, namely Lyapunov exponents, is difficult to employed. Therefore, various, alternative diagnostic methods have been proposed, as for example, 0-1 test [12,13] and the method of wandering trajectories [10,11].
This work investigates the opportunity to apply recurrence analysis that has not been utilized yet for SMA oscillators. This approach includes Recurrence Plots (RP) and Recurrence Quantification Analysis (RQA) that have been of interest to characterize the nonlinear dynamics of many dynamical systems [14], including solutions and bifurcation identification in mechanical systems [15,16] and smart material structures [17]. In this method, the periodicity of the underlying dynamics can be investigated by means of recurrences that are calculated for each point of the reconstructed trajectory. This method was originally developed by Eckmann et al. [18] and directly extended by Webber and Zbilut [18,19], Casdagli [20], and finally, by Marwan et al. [14]. The application of recurrence analysis can be applied for both short deterministic and noisy experimental data (see [15,21]).

Modelling of the shape memory alloy oscillator
An oscillator, composed by a mass m and a restoring force provided by a superelastic SMA device, is considered. The oscillator is subjected to a harmonic excitation represented by F(t) (Figure 2). The thermomechanical constitutive model, proposed in [22], has been adapted in [23] as a model for the restoring force of SMA devices, suitable for the use in the analysis of nonlinear dynamics of Shape Memory Oscillators (SMO) [10,11].
This model is employed to form the SMO analysed in this paper, considering a dimensionless form discussed in the sequence. Specifically, quantities with the dimensions of a temperature are normalized with respect to a reference value θref, greater than the transformation temperature Af, typical for each specific SMA. On the other hand, quantities with the dimensions of time, length and force are normalized, respectively, to the elastic fundamental period 1/ω of the oscillator and to the displacement xMs and force fMs at the beginning of the forward A  M transformation at the temperature θref. The SMO is subjected to a harmonic excitation Φ(t) = FcosΩt, and operates in a thermomechanical environment, where a convective rate of heat exchange Q(t) = h(θe-θ(t)) takes place. Here, F and Ω are the normalized excitation amplitude and frequency, θe the environment non-dimensional temperature and h the coefficient of convective heat transfer.
At each time t, the state of the SMO is described by displacement x(t), velocity v(t) and temperature θ(t) of the device, as well as the martensitic fraction ξ(t)  [0,1]. Moreover, to model the complex hysteretic behaviour of SMA related to internal subloops, the state of the device also depends on the value ξ0(t), that represents the martensitic fraction ξ, attained at the beginning of the last phase transformation, occurred before time t. As discussed in [8] Bernardini and Rega (2005), time evolution of the SMO is described by the following system of 5 differential equations: (1) The constitutive functions Λ and Z take different expressions, depending on the fulfillment of suitable transformation criteria that detect the occurrence of the phase transformations [8]. Subscripts F and R indicate the expressions relative to the upper and lower plateaus, respectively, that correspond to the forward (A  M) and reverse (M  A) transformations: where ΨF,R are constitutive functions that determine the shape of the upper/lower pseudoelastic plateaus and the other quantities are model parameters. Whereas various choices at different levels of refinement are possible, the following one is adopted: Here b is a parameter that determines the smoothness of the transition between the elastic branch and the transformation plateaus. The second law of thermodynamics expresses the non-negativity of the rate of energy dissipation Γ. In this framework Γ=Λ .
  Since Λ=ΛF when  >0 and Λ=ΛR when 0    , the second law requires: ΛF ≥ 0, ΛR ≤ 0. Such constraints impose limitations on the range of variations of model parameters. For example, while evaluating ΛF at the beginning of the upper plateau, the following constraint among λ, J, q2 is obtained: Considering (1) and (2), the system response depends on 7 model parameters, besides the above-mentioned b: q1, q2, q3, λ, L, h, J.
The first four parameters determine the shape of the outer pseudoelastic loop. In particular: q1 and q3, On the other hand, q2 determines the position of the isothermal lower plateau with respect to the upper one hence, determines the size of the hysteresis loop. Finally, λ directly influences the length of both plateaus. The remaining three parameters characterize the thermomechanical properties.
In particular, a parameter L determines the amount of heat which is produced or absorbed during mechanical loading, and the parameter h characterizes the rate at which the involved heat can flow out from the system to the environment by convection. Physically meaningful ranges of values may be identified with L  [0.0, 0.5] and h  [0.0 ,0.2]. The limit case L = 0 corresponds to shape memory materials which transformations produce negligible amounts of heat, while the limit case h = 0 models an adiabatic environment in which all the heat produced remains in the system.
The thermo-mechanical parameter J determines the slope of the linear dependence of the transformation forces on the temperature. A physically meaningful range of values may be identified with J  [1.0, 4.0].

Recurrence analysis
In early beginning, the concept of recurrences was applied to study the nonlinear dynamics and bifurcations in terms of stroboscopic phase diagrams, namely Poincare maps [24]. The method of the Recurrence Plot (RP), introduced by Eckmann et al. [17], was generalization on all the states along the trajectory [14,24]. RP is based on identification of the same states in the phase space of the dynamical system. The recurrence of the state xi at time i, at a different time j, xj ≈ xi, is counted within a two-dimensional squared matrix R with ones and zeros, i.e.
where N is the number of considered state states xi, along the embedding space trajectory, defined in a discrete sampling time ti=iΔt, ε is a threshold distance, || || indicates a norm (usually the L2 norm) and Θ (.) represents the Heaviside step function.
The graphical representation of such a matrix, where ones and zeros are represented by black and white dots respectively, is called recurrence plot (RP). In other words, information provided by the RP can be described as follows: 1: , , 1, , , 0 : where xi ≈ xj indicate the points from the neighbourhood of radius ε (defined according to the applied norm). In order to calculate the RP representing properties of the considered system, it is necessary to determine the embedding space and define values threshold ε. Embedding parameters include the time delay, τ, to introduce missing independent variables as vectors xi =[xi(1), xi(2), xi (3), …, xi(m) ], where xi(n) = xi-(n-1)τ with the appropriate embedding dimension, m, where n ≤ m [26].
To estimate the smallest sufficient embedding dimension m, the false nearest neighbours' algorithm is frequently used [27]. Appropriate time delay τ can be determined with the use of the auto-correlation or mutual information function [28,29]. In the original definition of the method, the neighbourhood is specified with the use of the L2 norm and its radius is chosen in such a way that it contains a fixed amount of states. For such a vicinity, the radius εi changes for each xi (i = 1, …, N). This neighbourhood is often denoted as fixed number of nearest neighbours (FAN). Other commonly used neighbourhood is that of a fixed radius, firstly used by Zbilut and Weber [30].
During research, carried out over the last few years [15,[31][32][33], it has been proven that the recurrence-based methods are particularly sensitive to changes in the dynamic behaviour of mechanical systems. In this paper, the methodology based on the recurrence plots (RP) method was applied to analyse regular and chaotic motions of the superelastic shape memory oscillator. In all computations, FAN definition of the neighbourhood was used [25]. The smallest sufficient embedding dimension (m = 3) was estimated with the application of the false nearest neighbours' algorithm [25]. Appropriate time delay (τ = 11) was determined by searching for the first minimum of the mutual information function [29], while the phase space of the considered systems has been reconstructed by delay embedding [30,34]. For the purposes of consistency, for both signals, the same value of threshold (ε = 0.09) was assumed. All the necessary computations were performed with the use of the CRP Toolbox for MATLAB [35].
Recurrence Quantification Analysis (RQA) [14,30] is a method of nonlinear data analysis which quantifies the number and duration of recurrences of a dynamical system, represented by its state space trajectory. Definitions of the most popular RQA measures, such as: The Recurrence Rate (RR), Determinism (DET), Laminarity (LAM), Averaged diagonal line length (L), Trapping Time (TT), Longest diagonal line (Lmax), Longest vertical line (Vmax), Divergence (DIV), Entropy (ENTR) can be found in [30,36]. Extended discussion, including recurrence times of the first and second orders as T1 and T2, can be found in Marwan et al. 2007. On the other hand, the additional topological measures, as transitivity or causality (TRAN) and clustering coefficient (CC), were given in [37].
The simplest measure of the RQA analysis is the recurrence rate (RR): defined as a ratio of all recurrence states (points) to all possible system states. Therefore, RR describes the probability that a state recurs in its εneighbourhood in phase space. Definition of another measure -DET -is based on the observation that processes with uncorrelated or weakly correlated, stochastic or chaotic behaviour, cause none or very short diagonals, while deterministic processes result in longer diagonals. Therefore, the ratio of recurrence points forming diagonal lines to all recurrence points is treated as a measure for determinism (predictability) of the system:  (8) where lmin often is fixed to the shortest line length lmin =2 and P(l) denotes a histogram of diagonal lines of the length l: where Nl is the total number of diagonal lines. Depending on the system character divergence of the dynamical system can be associated 1/Lmax or 1/Lmean, where Lmean is an average length of diagonal lines, can be used to distinguish the chaotic and periodic systems. Periodic systems would have long lines and chaotic broad spectrum of line lengths resulting with increase of divergence.
The ratio between the recurrence points forming the vertical structures and all the recurrence points is called laminarity (LAM). A value of LAM decreases if the RP consists of more separated recurrence points than vertical structures: Trapping time, TT, expresses the average length of vertical structures: It can be interpreted as the time during which the system stays in a particular state.
Transitivity, TRAN, is a network topological parameter defined on the recurrence matrix, R, elements: In general TRAN is a measure which refers to the binary relations between three given points i, j, and k, and it is sensitive to the dynamical system attractor geometry [16,38]. On the other hand, the maximal length of the vertical structures in the RP, denoted as Vmax: is a measure, which is the analogue to the standard RQA measure Lmax (here, Nv denotes the absolute number of vertical lines). In contrast to the RQA measures based on diagonal lines, the measures based on vertical lines can find chaos-chaos and chaos-order transitions.
Furthermore, the recurrence times of the first, T1, and of the second, T2, type define the average time distances between the recurrence structures. T1 includes the average horizontal space between subsequent structures, in terms of the white points, while T2 corresponds to the full periods including white and black points. Both parameters give an idea on system periodicity given in RP by the lengths between (basically diagonal) lines.

Results
As discussed in detail in [8], the SMO can undergo chaotic motions. For the purposes of the presented comparison, two specific trajectories are shown: one chaotic and the other -periodic. Both trajectories are obtained after numerical integration of the system for 500 periods by means of a standard 4 th order Runge-Kutta algorithm with 4000 steps per periods. The output is then sampled to get time series of 200,000 points. Model parameters are summarized in Table 1 corresponding to devices with a hysteresis loop of medium-high size and low hardening. The latent heat of transformation is low so that thermal phenomena are small, and the resulting conditions are close to the isothermal ones. The two trajectories differ for the features of the twoforcing excitation (Tab. 1). The simulation results using the above presented model are given in Figure 3. Time histories and corresponding spectra are plotted in Figures 3a, c, and 3b, d for periodic and chaotic signals, respectively.  Note that Fourier spectrum of the periodic response is given by a singular spectrum (Figure 3c), while the chaotic response is characterized by a continuous spectrum with additional singular peaks of higher amplitude. Occurrence of higher singular periods (with respect to excitation frequency Ω/2π) in the spectrum is a consequence of system nonlinearities (see Eq. 1).
The spectrum criterion has a qualitative character. A standard way to categorize such solution would be using Lyapunov exponent or 0-1 test [12]. In the sequence, the recurrence analysis is performed, showing that it can be used as an alternative way to estimate the character of solutions [14,15].
The recurrence analysis includes recurrence plots (RPs) and estimation of the recurrence quantifiers (Equations 7-15). The RPs are presented in Figure 4. Note, that periodic response is characterized by the collection equally spaced diagonal straight lines ( Figure  3a). On the other hand, the chaotic response is based on the shorter curved lines with evident variable distances between them. This basic qualitative difference included in the RPs is mirrored in the recurrence parameters (Table 3). Namely, values of particular RQA measures, estimated for both cases, represent quantitative information sufficient to decide about the character of system response. Table 3. RQA measures computed for the considered periodic and chaotic signal for RPs from Figure 4. Assumed parameters: m = 3,  = 11, ε = 0.09, FAN criterion. Calculations were done using the numerical package [35]. This occurs due to the chaotic solution being complemented by some periodic features. These features are related to the coexistence of the broad band of frequency and narrow peaks. This structure is mapped into RP (see Figure 4), where the line layout in Figure 4b follows, in some extend, the line structure of Figure 4b. However, it is possible to observe important deviations into the discontinuous and curved lines. This specific behaviour can show the intermittency on the system [16]. Different recurrence parameters support this possibility.

State
Firstly, we notice the difference in Lmean. As expected, chaotic response shows larger divergence (defined as 1/Lmean). Secondly, entropy ENTRL is 3 times bigger for chaotic system than regular one. Thirdly, Vmax is larger for chaotic system (4 times larger); this is due to the appearance of the curved lines characteristic for variable frequency responses. This result, together with the results of the trapping time, TT, in much larger for chaotic system support a scenario of system persistence around a specific state mediated by more periodic motion. Finally, the additional arguments are provided by the T2 and TRAN parameter. Obviously, it can be noticed that the ratio T1/T2 indicates on some smearing in response phase (local decelerations or accelerations in the system response). TRAN shows the information about the reproduction of binary related between black points on RP. It is obvious, that the complexity of the studied system responses is different.

Conclusions
We have studied nonlinear dynamics of a pseudoelastic shape memory alloy oscillator. The hysteretic behaviour of the SMA element is described by a phenomenological model with internal variables. Numerical simulations are performed showing periodic and chaotic responses. Diagnostic analysis of the kind of response is performed considering recurrence plots. Recurrence plots and recurrence quantification analyses are considered. Using various recurrence measures, a successful distinction of regular and chaotic solutions has been obtained using the following parameters: Lmean, ENTRL, TT, T2, Vmax. Results are compared with Fourier spectra presenting coherent results. In general, it is possible to conclude that the recurrence analysis is capable of distinguishing periodic and chaotic responses from short intervals (ten cycles of excitation -see Fig. 4).