Abstract
Complex systems are a broad concept that comprises many disciplines, including engineering systems. Regardless of their particular behavior, complex systems share similar behaviors, such as synchronization. This paper presents different techniques for determining the source of coupling when a set of oscillators synchronize. It is possible to identify the location and time variations of the coupling by applying a combination of analytical techniques, namely, the source of synchronization. For this purpose, the analysis of experimental data from a complex mechanical system is presented. The experiment consisted in placing a 24-bladed rotor under an airflow. The vibratory motion of the blades was recorded with accelerometers, and the resulting information was analyzed with four techniques: correlation coefficients, Kuramoto parameter, cross-correlation functions, and the recurrence plot. The measurements clearly show the existence of frequencies due to the foreground components and the internal interaction between them due to the background components (coupling).
1. Introduction
Complex systems are common in engineering, examples being communication networks, computer networks, and electronic devices, as well as mechanical systems such as manufacturing plants, aircrafts, and automobiles. Determining the dynamic behavior of these mechanical systems requires the use of advanced analytical techniques. As the number of elements increases, simple systems evolve into ones that are more complex. Concepts such as correlation dimension, Lyapunov exponents, fractal properties, and complex network analysis are some of ideas that are used to determine dynamic properties. The problem is to identify the source of synchronization from field data.
In addition to other behavioral characteristics, complex systems can present synchronization. In most cases, the system deploys two behaviors: one that is easily identified at a foreground level and another that remains at background level [1]. The foreground behavior is seen at frequency responses with high amplitudes and narrow bands, whereas the background behavior presents a wide band response and low amplitudes. The study of synchronization has been an exciting and interesting problem since Huygens first wrote about his experiments in 1665 [2]. Peña Ramirez et al. [3] analyzed Huygens’ experiment, modeling the system as a lumped mass model, and found that synchronization depends on the relationship between the masses, the stiffness of the connecting beam, and the amount of damping in the system.
Keeffe et al. [4] presented a model for the analysis of pulse coupled oscillators. As they described it: in some systems synchronization starts at a certain location and expands forming clusters. They suggested that clusters formed, evolved, and collapsed. For this purpose, they developed a model that provides an idea of the physical phenomena and how clusters evolved into synchronization. They described different systems that behave as pulse coupled oscillators. They presented two types of oscillators, systems with local coupling and systems with global coupling. There is still a question how transient dynamics leads to synchronization [4]. Ulrichs et al. [5] presented the analysis of metronomes as nonlinear periodic oscillators. The metronomes behave asymptotically and they show higher dimensional attractors. Metronomes have the same behavior as Huygens clocks, and their synchronization can be modeled using Kuramoto’s parameter. They found small traveling waves through the system support that lock the phase difference among all metronomes and they model the connections in terms of a Van der Pol oscillator. Aragonès and Guasch [6] presented a theoretical analysis of path propagation in a set of systems with similar structures. Woodhouse [7] and Langley et al. [8, 9] introduced statistical energy analysis. They modeled the system as a set of springs and represented it as a block diagram. The diagram represented subsystems as blocks and power flows as connections among subsystems. For the net analysis, they assigned weight factors relating a coupling loss factor to a total loss factor. Moreover, they dealt with the net analysis including a probabilistic density function and computed the weight factors as a statistical ratio. They defined the length of a particular link as a random variable, and they defined this function as a mean function plus a function of the standard deviation. A similar approach was presented by Llerena-Aguilar et al. [10]. They reported a hybrid method for the analysis of wave propagation and synchronization of a Wireless Acoustic Sensor (microphones) Networks. The problem they solved is the inverse of what we are trying to understand. In fact, they studied the desynchronization. The phenomenon is caused by two factors, namely, the clock phase offset of each sensor and the propagation delay due to the distance between them. These factors are critical in the selection of the processing algorithm. The classical solution was the implementation of clock synchronization protocols. The phenomenon is simulated assuming that the source signals are received by an array of sensors. Then, the signals are combined adding a Gaussian noise and convoluting the signal with an impulse function for each sensor. They discretized the signal using the Short Time Fourier Transform and they aligned the mixture of signals using the cross-correlation and the Log Energy Cross-Correlation functions. Ray et al. [11] investigated the Noise Induced Synchronization in two identical time delay systems. They compared numerical results with experimental data and they analyzed the effect of white, green, and red noises. Messina and Vittal [12] reported a method for analyzing dynamic patterns from measurements. They applied the Hilbert transform and proper orthogonal analysis technique. This technique is based on the relation of the energy and phase embedded into the mode function. Rosenblum and Pikovsky [13] presented a method for detecting weak couplings between two oscillators. If two oscillators are weakly coupled, they assume that the two signals contain enough data from the weak interaction. The relationship in a driven-response from noisy data could be unreliable and the interpretation is difficult. They discussed the effect of noise as follows: “in case of perfect synchrony, we are not able to separate the effect of interaction from the internal dynamics of autonomous systems. In order to obtain the information on the direction of coupling we need to observe deviations from the synchrony, either due to noise or due to onset of quasiperiodic dynamics outside the synchronization region.” They propose a method for determining the direction of the coupling and they quantified the cross-dependencies of the phase dynamics. They illustrate their method with two coupled Van der Pol oscillators.
Friedrich et al. [14] presented an approach for the study of complex systems using stochastic methods. Complex systems can be traced back to rather simple laws because of the existence of self-organized processes. While the analysis of nonlinear dynamics is traditionally conducted using time series, complex systems show stochastic data in time and length scales. It is clear that the synchronization of complex systems is a function of time and space. In Friedrich’s work, the dynamic model is represented as a Langevin equation with white noise. Since the dynamic model can be represented as a deterministic set of nonlinear differential equations, they added a stochastic term with white noise. They considered a Gaussian distribution for time propagation as a Markov process. A model for self-similarity was included, meaning that at a certain time two elements should have the same statistics. In a previous work, Ghasemi et al. [15] presented the steps below for developing the stochastic equation. First, review if the data follow a Markov chain. Then, estimate the Markov time scale, which is the minimum time interval for a Markov process. For this analysis, it is necessary to use a joint probability density function. This can be simplified with the Chapman-Kolmogorov equation, which is the probability of having two similar distributions in a time interval. Then, they derived the stochastic equation and a regeneration of the stochastic process. They applied the method to heart beat signals. In complementary analyses, Siefert et al. [16] and Shirazi et al. [17] classified complex systems as networks with fluctuating stochastic processes. For this problem, they y represented spatial networks of complex dynamics in phase space. They proposed a method for mapping a stochastic process onto a complex network. As for recurrence plots, they determined the Markov time scale, and they based their analysis on the Markov process.
Marwan et al. [18] presented a review of recurrence plots applied to dynamic systems. Recurrence is one of the basic properties of dynamic systems and the recurrence plots can help visualize dynamic system behavior. They described the use of recurrence plots for identifying synchronization in the sense that two systems synchronize if their recurrence in phase space is linearly dependent. In another work, Donner et al. [19, 20] studied the dynamics of complex systems using the recurrence concept and network theory. They demonstrated that exploiting the duality of the recurrence matrix and the adjacent matrix of a complex network improved the determination of quantitative characteristics of recurrence networks. They divided them into three types of recurrence networks, that is, local, intermediate, and global. All their results were obtained from simulations.
Other researchers have applied different techniques to identify synchronization. Huo et al. [21] analyzed the synchronization of a couple of calcium oscillators. They proposed a mathematical model based on a set of nonlinear first-order differential equations. They found synchronization through the analysis of phase diagrams. Kim and Lim [22] worked on burst and spike synchronization of neurons. Their interest was aroused by the analysis of electroencephalogram (EEG) data. Researchers believed that synchronization of neuron activities could help in understanding sensory and cognitive processes. From an engineering point of view, it can be understood as periods of high frequency impulse generation and periods of no impulse generation. Xu et al. [23] and Rosário et al. [24] described motif synchronization. They studied temporal and spatial synchronization from EEGs and applied time-varying graphs, static network, and motif synchronization. Motif synchronization is a sequence of small patterns in a particular order of occurrence. With three first-order differential equations, they simulated the propagation of an external stimulus in time and space.
Childs and Strogatz [25] analyzed a periodically forced Kuramoto model. They defined a rich dynamic system when it has a strong interconnection among randomness, force, and coupling. Randomness is associated with variations in the natural frequencies; this concept is important because it determines the desynchronization of the oscillators. The coupling is related to the oscillator’s phase and the force to the external excitation. Thus, the system synchronizes and desynchronizes depending on the relative magnitude of these three effects. Another important concept described by Childs and Strogatz is that a subset of the oscillators is “phase-locked” to the force drive, whereas the rest of the oscillators are on the tails of the frequency distribution outside of synchronization. Their method is applied to a set of Van der Pol oscillators and they proposed that the oscillators are sinusoidal coupling; in our case it is assumed that the coupling has a very low frequency compared to the natural frequency and the amplitude is low.
The objective of the present paper is to determine, using experimental data, if there is a way to identify the background effects that enable synchronization in complex systems and coupling. It is assumed that a coupling source exists and is the cause that synchronizes all the elements. The difficulty of identifying this coupling source rests on the fact that its dynamic response is of low frequency and low amplitude. First, different techniques that describe synchronization are applied; then, data evolution is analyzed in order to determine if there is a very slow dynamic behavior. The results show that synchronization is easily identified, but the source of synchronization is not; the reason is its low energy content and its larger fundamental frequency. Direct FFT, or other traditional dynamic analysis techniques, is unable to process this type of dynamics. The background modes are, however, responsible for internal interaction, communication, and energy transfer between the natural frequencies of the foreground components. The background vibrations are transmitted through somewhat undefined neighboring linkages at very low frequencies, which may be orders of magnitude lower than the natural frequencies of the foreground element.
2. Experimental Setup
The experimental setup consisted of a stationary rotor with 24 blades. In order to reduce the number of variables, the blades were built with a rectangular cross section, as shown in Figure 1, and they were mounted with the thin side of the blade facing the front. At the tip of each blade, a MEMS accelerometer was mounted; data from 22 accelerometers were recorded. These accelerometers were double-axis Analog Devices ADXL321 type with a sensitivity of 10 g. Their masses were small enough so they had a negligible contribution to the natural frequency of the blades. A laser Doppler vibrometer, similar to that used by Oberholster and Heyns [26], was used to calibrate the individual accelerometers that were mounted on each blade. Figure 2 shows the experimental setup and the laser vibrometer used for calibration. The accelerometers were connected to a data acquisition system, and the data were recorded simultaneously.
The rotor was tested under two types of loads, namely, an impact and an airflow test [27]. For the impact test, a calibrated hammer was used to hit the center of the rotor and all the blades vibrated simultaneously. Data were collected for 1.6 s at a sampling rate of 5000 Hz. The airflow test consisted of blowing air through a duct producing an airflow perpendicular to the blades. The flow velocity was kept constant and accelerometer data were recorded continuously over a time interval of 2 s at a sampling rate of 1 kHz. There were ten recordings intervals, called events (Table 1). Event A started when the airflow was switched on.
3. Data Analysis
The natural frequencies of each blade were determined with an impact test. Figure 3 shows the Fourier spectrum of blade 1. It is clear that there are two dominant peaks related to the foreground components, 120.9 Hz corresponding to the natural frequency of the blade and 59.9 Hz to the rotor structure.
The blades were found to have almost the same natural frequency between 121 and 122 Hz. Since the present analysis is focused on determining the state of synchronization, for the airflow experiments each blade response was filtered around 120 Hz. Figure 4 shows the frequency spectrum of the airflow test corresponding to event A. In this way, only the natural dynamic response of the blades was included in the analysis. Filtering is an important step since other vibration sources hide synchronization, and the lock to the driven force is minimized [25]. It was found that the analysis techniques are very sensitive to external sources of vibration. The analysis presented in this paper is centered on the difference between each foreground element; therefore, if there is a strong source of vibration, every accelerometer will measure the same value and the contribution of the interaction between the background components will disappear.
4. Correlation Coefficient
The correlation coefficient is a method that calculates the angle between two vectors. This coefficient calculates the inner product of the data vectors and measures the linear dependency of two signals. Although it differs from the phase angle, it gives a good estimation of the synchronization of two signals. Applying the Pearson correlation coefficient to the 22 data vectors a symmetric 22 by 22 matrix is obtained. Unfortunately, since this parameter measures the overall phase angle, it cannot identify the internal interaction that brings the elements into synchronization.
The correlation coefficient for the impact test is displayed in Figure 5. This figure presents a pattern that is intrinsic to the actual structure. Ideally, this figure must have the same values except at the diagonal, but the differences correspond to deviations in the construction of each blade. This figure is the characteristic signature of the structure.
Every map showed a similar pattern to the impact test. This pattern is associated with the rotor structure; it is a signature of the system and appears regardless of the analysis. As a result, the correlations between blade 5 and blades 13 and 16 are greater than 0.5 in every map. Similarly, blade 10 has a strong correlation with blade 20. The rest of the blades have oscillatory values for the entire test. From these maps, it is difficult to define a function that describes these oscillations, since they vary in time and position within the matrix. Therefore, the magnitude of each matrix was evaluated. Figure 6 shows the variations of the accumulated value as a function of time:
It can be seen that after event C there is little variation.
The correlation factor is unable to identify synchronization, and to overcome this the data were analyzed as indicated next.
5. Kuramoto Parameter
The Kuramoto model describes analytically the phase transition of a set of oscillators into a mutual synchronization and determines a threshold at which the oscillators are coupled [25]. The Kuramoto parameter measures the instantaneous phase angle of several signals. It is based on the analysis of the synchronization of a set of oscillators moving at a similar frequency but with a variable phase. This parameter was used to determine the amount of synchronization among the blades. The instantaneous phase angle was determined using the Hilbert transform. This is given by where and , both, are real-valued time domain signals.
From the Hilbert transform, a complex vector is defined as
Thus, the instantaneous phase angle is calculated as
After calculating the instantaneous phase among the blades, the Kuramoto parameter was determined to be where , , represent the instantaneous phase angle.
When the Kuramoto parameter is equal or close to 1, all the oscillators move together. If it is equal or close to zero, then the oscillators move independently and no synchronization occurs. Since the actual data contain several frequencies, the Kuramoto parameter will evaluate the contribution of all frequencies, including the excitation frequency. In general, this frequency will dominate over the rest; thus, it was necessary to filter the response data around the natural frequency. With this procedure, the Kuramoto parameter was calculated at each event. To have an estimation of the degree of synchronization, the histogram of the Kuramoto parameter was calculated, and it was assumed that the blades are synchronized if the majority of the computed values are above 0.7. The Kuramoto parameter is a good indication of synchronization. Nevertheless, it is unable to identify the source of synchronization.
The result of analyzing the impact test is shown in Figure 7.
In every event, most of the time the Kuramoto parameter has a value around 0.7, indicating that the blades vibrate synchronously. Nevertheless, it is difficult to identify the source of synchronization from these plots.
6. Cross-Correlation Function
This technique of analysis consists of determining the cross-correlation function between the blades. This function is calculated between the acceleration data of one blade with respect to the other 21 blades. It estimates the statistical correlation between two signals at two different points in space or time. It is a useful indicator of dependencies or similarities in time of two signals. In addition, it can be used to interpolate for values at instants for which there are no observations. With this basis, the cross-correlation function determines similarities among the data and disregards those instances when the blades vibrate out of synchronization.
Mathematically, the cross-correlation function is defined as where is defined as the lag.
Based on the elements of a rich dynamics, the internal interaction between the background components occurs in time and space; the cross-correlation function can be used to produce maps in time and space; time is represented by the lag order, being a normalized parameter that is proportional to time, and the blade number represents space. At each event interval, there are 22 maps. The maps are constructed by where is the reference blade and are the remaining 21 blades.
In this way, discrepancies among similar signals can be identified in time and space. Since the function is symmetric, only half of it is used in the analysis. With the construction of these maps, it is possible to identify three independent variables: a short time (lag order), a long time (defined in Table 1 as event), and the space (relative position between each blade and blade 1). Based on these three variables, contour plots were produced for each measurement data. Figure 8 shows the contour map of the impact test with respect to blade 1. The horizontal axis is the blade number in (7), the vertical axis is the lag order (), and the contour plot represents the magnitude of C.
As in the correlation factor, there are three dominant regions or clusters (vertical strips) at blades 6, 14, and 22. From the contour plots, it is easily noticed that the correlation amplitude increases linearly with respect to lag. Nevertheless, there are variations at certain blades.
After analyzing the impact test, the cross-correlation functions of the ten airflow measurements were calculated. These calculations produced 10 × 22 maps where it is possible to see variations in time and space. As an illustration, the ten maps corresponding to blade 1 () are presented in Figure 9.
The maps show three different variations, the maximum correlation value, the density of contour plots, and the location of the maxima. This suggests that synchronization varies in time and space. Determining these variations is not straightforward, and for that reason the four different analyses were conducted. The first consisted of evaluating the variation of the maximum correlation value. Figure 10 shows the blade number where the maximum correlation value occurred, each curve corresponding to an event. Most of the maximum values occur at the same blade location and they evolve following the same pattern, except for blade 10. Figure 11 shows the evolution of the maximum correlation of blades 6, 10, and 19. The maximum has an incremental tendency as the test time increases, whereas the maximum correlation of blade 10 decreases. After comparing these figures, it seems that the internal interaction of the background components starts at blade 10 and propagates to the adjacent elements.
7. Recurrence Plots
Another technique of analysis that can identify synchronization is the recurrence plot. Recurrence is a concept related to “memory,” if the behavior of a system can be predicted; it means that the system returns to a former state. Poincaré proposed the original concept, and Eckmann introduced the concept of recurrence plots. In general, from Hamilton’s principle, the trajectory along the phase space depends on the velocity and the displacement. The phase space of a simple harmonic has the shape of a helix along the time axis direction. According to Marwan et al. [18], the recurrence plot determines the evolution of the system and finds the number of times that a signal passes through the same point in the phase diagram. If two signals are synchronized, their phase diagrams would be identical. The recurrence plot is a graphical representation of a matrix formed with all the recurrence values. A pure harmonic system shows evenly spaced diagonals. However, for nonharmonic systems, the shorter the diagonals the more chaotic the system. Marwan et al. probed the theory using the Poincaré recurrence theorem. Nevertheless, the theorem is insufficient for predicting the time for a system to return, that is, the recurrence. The determination of the time interval has similar properties to the Markov process, and it is related to the Lyapunov exponents and entropies. What they presented is a way of analyzing the length of the diagonals. One of the challenges of this method is the determination of the appropriate time interval, which is the difference between an orbit and the following orbit. The phase trajectory is represented as an evolving path; thus, the time evolution is the derivative of the trajectory with respect to time. In their work, the authors proposed a method to reconstruct the phase space from discrete data with a time delay factor. In practice, this is feasible only when the system response has a single dominant frequency; otherwise, it is necessary to use other techniques to reconstruct the third coordinate. To calculate the matrix, they recommended the use of the maximum norm. Since the recurrence plot is determined from the difference between a point and the subsequent point on the phase space, they proposed the use of a tolerance error for computing this difference.
In this work, it was found that the construction of the phase diagram from measurement data is cumbersome. Since the data are the acceleration of each blade, they have to be integrated numerically in order to obtain speed and position.
In order to analyze the synchronization between two signals, the data from the two vectors are analyzed at the same time span, the two phase diagrams are superimposed, and the distance between every point is calculated. In this way, the recurrence plot is constructed using the equation
Vectors and are normalized using the maximum value for each case. Ideally, the value of should be zero, but in reality it is impossible. The criterion for determining synchronization is that , where is a tolerance value.
As in previous analyses, the first plot corresponds to the impact test. Then, each event is analyzed. Figure 12 shows the recurrence plot of the mean and the standard deviation. The mean value is calculated taking each vector as a reference, starting with blade 1 and finishing with blade 22. In this way, a 22 × 22 matrix is obtained: where is the reference vector and is the corresponding vector, is the instantaneous value of vector , and is the instantaneous value of vector .
(a)
(b)
It is interesting to notice that the structure of the plots is different. With this analysis, the impact test produces a plot that has a different structure than the other data. In the entire airflow test, the dominant values correspond to blades 6, 14, and 21. Data create clusters around these blades, and the cluster pattern appears in every map, except during the impact test. Nevertheless, the pattern shows variation in every event. These variations are associated with the slow dynamic behavior that is a manifestation of the internal interaction between the background components.
In order to determine the slow dynamic behavior, the sum of individual values is calculated as
This number reflects the amount of synchronization at each event. Ideally, it should be zero, but since there are variations in the individual phase diagrams, it reflects the overall dynamic behavior. Figure 13 shows the variations of as a function of time.
This figure shows a slow cycle with a period of about 36 to 40 min. Therefore, the internal interaction between background components has a frequency of Hz. With this analysis, it was possible to estimate the time variation due to the internal interaction. This depends on time and location; the time variations were identified with the recurrence plot and the location with the cross-correlation function.
8. Conclusions
In this paper, the study of the synchronization of a rotor with 24 blades is presented. 22 acceleration vectors were registered simultaneously, and the rotor was subjected to two types of tests: an impact test and an airflow test. Data were recorded at ten different events, separated by long intervals of time.
The measurement and analysis of field data present different challenges; the length of the data is limited by the data acquisition system, and therefore it is impossible to record continuously the acceleration of each blade over a long period of time. The other challenge is the waveform of the internal interaction; in general, the internal interactions have a very low energy content and very low frequency. The accelerometers are unable to identify these waveforms; thus there is a need for other techniques of analysis that are able to identify this behavior.
Synchronization can be easily identified with the correlation coefficient or with the Kuramoto parameter. It was demonstrated that these techniques show synchronization among the foreground components, but they were unable to establish the internal interaction between the background components, neither in time nor in space. In order to overcome this limitation, the data were analyzed with the cross-correlation function and the recurrence plot.
With the cross-correlation function, the data were processed in a way that variations in time and space were represented as contour maps. These contour maps were constructed with the local time variable, the data source (the blade number), and the long waveform (the ten events). With the cross-correlation function, the source of the internal interaction was identified. The criterion for identifying the source was the discrepancy of the behavior of blades 6 and 19 with respect to the rest of the data. These blades showed greater correlations at early events than the rest of the blades that evolved following a similar pattern.
The recurrence plot was the technique that allowed the analysis of the time evolution. This was constructed from the phase diagrams of the blades. Taking the difference of the phase diagrams among the blades it was possible to determine the amount of synchronization. With these differences, the evolution of the mean value during the ten events was evaluated. The evolution in time showed a very slow dynamic pattern with a very low frequency. This frequency is associated with the time variation due to internal interaction among background components.
It has been shown that it is possible to identify the dynamics of internal interaction in complex systems. This interaction has a very low frequency and energy content; therefore, special analysis should be done. The correlation enables identifying the structure of the system and the evolution of synchronization; Kuramoto’s parameter estimates the amount of synchronization. The cross-correlation function allows identifying the location of the source, whereas the recurrence plot identifies the time variations. These analysis techniques provide some insight into the coupling source if they are applied together.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.