Synchronization in arrays of coupled self-induced friction oscillators

We investigate synchronization phenomena in systems of self-induced dry friction oscillators with kinematic excitation coupled by linear springs. Friction force is modelled according to exponential model. Initially, a single degree of freedom mass-spring system on a moving belt is considered to check the type of motion of the system (periodic, non-periodic). Then the system is coupled in chain of identical oscillators starting from two, up to four oscillators. A reference probe of two coupled oscillators is applied in order to detect synchronization thresholds for both periodic and non-periodic motion of the system. The master stability function is applied to predict the synchronization thresholds for longer chains of oscillators basing on two oscillator probe. It is shown that synchronization is possible both for three and four coupled oscillators under certain circumstances. Our results confirmed that this technique can be also applied for the systems with discontinuities.


Introduction
The synchronization has been widely studied in various fields of natural sciences. The classical experiment by Huygens back in the 17th century, who described the synchronization of pendulum clocks coupled by a common support, was the first scientific report on that phenomena [1]. Pikovsky et al. define it as an "adjustment of rhythms of oscillating objects due to their weak interaction" [2]. Examples of synchronization has been found in many fields of science [3][4][5][6][7][8][9].
For the purposes of this paper, we are focused on complete synchronization (CS). The CS can happen in case of identical systems when two state trajectories x(t) and y(t) converge to the same value and continue in such a relation further in time [10] or in other words fulfilling condition: lim t→∞ x(t) − y(t) = 0. Master stability function (MSF) is a useful tool to determine synchronization thresholds for different coupling configuration of identical oscillators [11]. According to classical method one needs Lyapunov exponents, as well as eigenvalues of connectivity matrix to obtain the MSF. In case when it is difficult to evaluate the MSF using classical algorithm, a technique which relies on above mentioned, simple definition of CS, called three oscillators universal probe can be applied [12]. For mechanical oscillators, a three oscillators probe can be reduced to only two oscillators probe [13], since between such systems only mutual and symmetrical linkage, also called real coupling [12], can be realized. Wu and Chua formulated conjecture [13] where they have stated relation between synchronization threshold and largest non-zero eigenvalues of connectivity matrix for networks of two different lengths of similarly coupled systems. Having two networks of length m and n with respective synchronization thresholds σ m , σ n , and the largest non-zero eigenvalues γ 1(m) , γ 1(n) of connectivity matrices G m and G n one obtains: Some development of this approach allows us to explain an occurrence of discontinuous synchronous intervals in coupling coefficient space -the phenomenon called ragged synchronizability [14]. Friction is one of the basics phenomena in mechanics and is responsible for resistance of contacting surfaces to motion. It can be observed in many practical applications like: wheels, musical instruments, stick-slip oscillations in machining, robot joints, hard drives, breaks. Researchers have developed different friction models, starting from the basic Coulomb model [15] to more advanced models [16][17][18][19][20][21][22], which take into account various properties of friction. Overviews of friction models are presented in [23][24][25]. The choice of the model depends on the characteristics of the investigated system. There is no need to apply sophisticated approaches, when the relative velocity between contacting surface is constant. However, in case of a rich dynamics of the relative velocity between contacting surfaces, non linear models of friction are preferable. Many friction models include discontinuity caused by the signum function, as the direction of friction is opposite to relative velocity. This can be problematic for numerical routines. Researchers have found many approaches to overcome this obstacle and we choose to approximate signum with inverse trigonometrical function arctangent (more details in Sect. 2).
We investigate the synchronization properties of self-induced dry friction oscillators coupled in arrays of different lengths. An example of application of the MSF method to the case of such strongly non-smooth dynamical systems is demonstrated. The MSF technique requires some indicators of the synchronization tendency. In this context, Lyapunov exponents seem to be the best natural tool. However, calculation of these exponents in the presence of discontinuity is not straightforward. Hence, we take advantage of the property given by Eq. (1) in our research and use results for two oscillators as a reference for evaluating synchronization thresholds of longer chains of coupled oscillators. The two oscillators probe significantly simplifies the synchronization analysis, as it allows to avoid the calculation of Lyapunov exponents.
The paper is organized as follows. The mathematical model used in the research is explained in Sect. 2. Results of numerical study are presented in Sect. 3. Finally, we draw conclusions in Sect. 4.

Model
Consider a simple, one degree of freedom (1DoF), forced stick-slip friction oscillator presented in Fig. 1(a). The system consist of mass m on the conveyor belt moving with constant velocity v b . The equations of motion for that particular configuration can be formulated as follows: where k is the stiffness constant, U amplitude of kinematic excitation, Ω angular excitation frequency, F N normal load force (note that the weight of mass m is already included in F N ), v r relative velocity between the contact surfaces and v b velocity Overdots stand for derivatives with respect to nondimensional time τ .
In order to model dry friction exponential friction model is applied, which is given by Eq. (4) and plotted in Fig. 1 The shape of the friction model function includes the Stribeck effect [26,27]. The signum function used in Eq. (4) contains discontinuity when relative velocity is changing its sign, which is not particularly comfortable for numerical integration software and can be solved by using either switch model, which detects whether the oscillator is in stick or slip phase or by using continuous approximation of signum function Eq. (5) where κ 1. In this paper approximation method is chosen with coefficient κ = 10 5 . Having derived equation for single oscillators, let us consider an array of n friction oscillators, which is presented in Fig. 2. All items are coupled by a spring of stiffness σk ,where σ defines the overall dimensionless strength of the coupling. The dimensional equations for that system are omitted and non-dimensional equations in matrix form are given bellow: 2672 The European Physical Journal Special Topics  with the corresponding connectivity matrix: defining the topology of connections between the oscillators. Coupling coefficient σ defines the strength of the coupling. Synchronization for n coupled oscillators can be measured by the synchronization error e n , which is defined as: Synchronization error describes trajectory separation of oscillators in phase space. For the CS state, the synchronization error is equal to zero. It is important to remember that the springs, which connect oscillators with the nearest neighbour are not the only coupling factors. All oscillators are subjected to the same excitation amplitude and frequency, conveyor belt velocity and have the same friction properties (i.e. model, coefficients), which should ease the process of synchronization.

Results
In this section, results obtained by the numerical integration of the equations of motion derived in Sect. 2 are presented and discussed. We focus on the synchronization properties of the system by using the MSF and average synchronization error as the research tools. The system is integrated by authors' own program written in C++ using Boost Odeint library with adaptive step size Cash-Karp method with tolerance levels at: e abs = 10 −15 , e rel = 10 −12 . The section is divided into several subsections, depending on the length of respective array. The values of parameters used in the investigations are given in Table 1. Should we use different parameters, an appropriate note is placed either in the caption or in the legend of each figure.

Single oscillator
Let us begin with the results for single friction oscillator (see Fig. 1). Investigation of the single oscillator is performed in order to validate the system and get information about the type of motion. Bifurcation diagram for assumed set of parameters is presented in the Fig. 3, with kinematic excitation frequency ω as a bifurcation parameter. Poincare cross-section for the bifurcation diagram are performed each time the velocity of the oscillator changes it sign from "+" to "−". The investigated case illustrates a cascade of period doubling bifurcations, as a route to chaos scenario. One can observe a fractal pattern, and periodic windows in non-periodic region ω ∈ [1.3, 1.9].

Two oscillators probe
For two coupled items the only non-zero eigenvalue of the connectivity matrix is γ 1(2) = −2. This is an important case, since it serves as a reference probe for the synchronization analysis for longer chains of oscillators. In general, we can define the MSF as a surface of the largest transversal Lyapunov exponent (TLE) over the complex numbers plane (α, β) representing arbitrary value of γ, where α = σRe (γ) and β = σIm (γ). As we mentioned in Sect. 1, in mechanical oscillators one has mutual and symmetrical interaction yielding to real coupling, hence the MSF can be represented as a function of the real number α, such that If all the eigenmodes corresponding to the discrete spectrum of eigenvalues σγ i can be found in the ranges of negative TLE, then the synchronous state is stable. According to an idea of two oscillators probe these synchronous ranges of α are approximated (via relation (10)) by corresponding σ-intervals of zero synchronization error. We are particularly interested in average synchronization error vs coupling coefficients diagrams, which are presented in Fig. 4, followed by bifurcation diagrams to investigate the correlation between synchronization and type of motion. Complete synchronization occurs when synchronization error is equal to zero. For lower excitation frequency ω

Three coupled oscillators
For the case of three coupled oscillators MSF is used to predict regions of possible synchronization. The connectivity matrix G 3 with two non-zero eigenvalues γ 1(3) = −1, γ 2(3) = −3 is as follows: Let us consider the system with excitation frequency ω = 1.15, when the system is in periodic regime. In Fig. 5 MSF e II (α) is projected via eigenvalues of connectivity matrix G 3 on bifurcation diagram of average synchronization error for 3 coupled oscillators for ω = 1.15. For that particular frequency of excitation the system behaves periodically. Synchronization error e II (σ) is scaled by absolute value of largest nonzero eigenvalue of connectivity matrix G 2 (i.e. γ 1(2) = 2) in order to obtain e II (α) (1) mentioned in Introduction. Thus, knowing beforehand the synchronization thresholds for two oscillators probe, we are able to predict synchronization thresholds for longer arrays of oscillators. When the angular frequency of excitation ω is increased to ω = 1.55 the system behaves non-periodically (see Figs. 4(c), 4(f)). Basing on the two oscillators probe for ω = 1.55 a similar procedure is performed to plot Fig. 6. Synchronization for three oscillators probe occurs for narrow interval of sigma with threshold values σ 1 = 0.176, σ 2 = 0.196, as in this interval α for both eigenvalues of connectivity matrix are in the region when synchronization error e II (α) = 0.

Four coupled oscillators
In case of four coupled oscillators one obtains three non-zero eigenvalues of the connectivity matrix  which are as follows: γ 1(4) = √ 2 − 2, γ 2(4) = −2, γ 3(4) = −2 − √ 2. The procedure for drawing MSF projection is the same as in previous subsection, the only difference is we have one additional eigenvalue. Examples of synchronization in periodic and non-periodic regimes are presented bellow. In Fig. 7 (ω = 1.15) the system finds it easy to synchronize for weak coupling σ < σ 1 , where σ 1 = 0.136, as all three eigenvalues meet in the synchronization α-region for two oscillator test probe. When the system behaves in non-periodic manner for ω = 1.4, the synchronization is hard to achieve and occurs only in limited region of σ (see Fig. 8). Since the slope of α for eigenvalue γ 1 is so low, the window when synchronization is possible is very narrow (i.e. σ ∈ [0.153, 0.177]).

Conclusions
We investigated the synchronization properties of coupled self-induced dry friction oscillators. Systems of two, three and four oscillators coupled by linear springs were tested to determine the synchronization thresholds using the MSF approach, basing on synchronization error of two oscillators probe. Our results confirmed that this technique can be successfully applied for the networks of strongly non-smooth dynamical systems. The detailed numerical analysis reveals that complete synchronization is possible for each of the tested network. Synchronization is more likely to occur in case of periodic motion. For the non-periodic motion synchronization occurs for narrow intervals of coupling coefficient σ. We show examples proving that synchronization threshold follow the MSF for two oscillator probe, although the system is also coupled by excitation and friction force. However, for longer array of oscillators (n ≥ 5), the complete synchronization regions disappear. For larger n, increasing ratio γ n−1 /γ 1 causes that not all values σγ i are located within synchronous α-range and then synchronization of all items in the array is impossible. On the other hand, the instability of synchronization (even for short array) may be due to possible high number of attractors coexisting in the phase space of the systems with discontinuities.