Nonlinear Dynamics of a Duffing-Like Negative Stiffness Oscillator: Modeling and Experimental Characterization

In this paper, a negative stiffness oscillator is modelled and tested to exploit its nonlinear dynamical characteristics. (e oscillator is part of a device designed to improve the current collection quality in railway overhead contact lines, and it acts like an asymmetric double-well Duffing system. (us, it exhibits two stable equilibrium positions plus an unstable one, and the oscillations can either be bounded around one stable point (small oscillations) or include all the three positions (large oscillations). Depending on the input amplitude, the oscillator can exhibit linear and nonlinear dynamics and chaotic motion as well. Furthermore, its design is asymmetrical, and this plays a key role in its dynamic response, as the two natural frequencies associated with the two stable positions differ from each other. (e first purpose of this study is to understand the dynamical behavior of the system in the case of linear and nonlinear oscillations around the two stable points and in the case of large oscillations associated with a chaotic motion. To accomplish this task, the device is mounted on a shaking table and it is driven with several levels of excitations and with both harmonic and random inputs. Finally, the nonlinear coefficients associated with the nonlinearities of the system are identified from the measured data.


Introduction
Devices and materials exhibiting a negative stiffness phase are often used as vibration isolators due to their amplified damping properties [1,2]. In particular, in the case of engineering structures, such devices are usually designed adopting discrete macroscopic elements, such as postbuckled beams, plates, shells, and precompressed springs, arranged in appropriate geometrical configurations. Examples can be found in automotive suspensions [3,4] or seismic isolation [5,6].
While the practical applications of negative stiffness systems are relatively recent, the theoretical studies about their dynamical behavior have a much longer story. is is because of the wide kind of motions they can exhibit, ranging from linear to highly nonlinear and chaotic [7]. In particular, when the negative stiffness effect is coupled to a nonlinear polynomial stiffness contribution, the so-called double-well Duffing oscillator is retrieved.
is oscillator exhibits two stable equilibrium positions plus an unstable one, and the oscillations can either be bounded around one stable point (in-well or intrawell, small oscillations) or include all the three positions (cross-well or infrawell, large oscillations). In both cases, periodic oscillations can evolve to steady in-well or cross-well chaotic motions under external excitation [8,9]. e occurrence of irregular motion, consisting of random-like crossings from oscillations around the two stable equilibrium positions, was first observed experimentally in 1971 [10] and the motion was called "snap-through". is behavior is associated with an unstable phase generated by the negative stiffness effect; thus, the stability analysis plays a crucial role in the design process of these kinds of systems [11]. Many studies have been conducted starting from the 1970s to analyze the dynamics of this oscillator and the effects of the unstable paths, including phenomena such as bifurcations, subharmonics, period doublings, and chaos. A comprehensive literature review can be found in [9].
In this work, a negative stiffness oscillator is studied and tested. e oscillator is part of a device designed to improve the current collection quality in railway overhead contact lines, attempting to alter their damping distribution and reducing the wave propagation [12].
A two-fold objective is pursued: first, the dynamical properties of the oscillator are analyzed, replicating experimentally the possible kind of motions it can exhibit (in-well, cross-well, and chaotic); second, nonlinear system identification is performed to extract the model parameters directly from the measurements. e latter in particular seems to be quite a challenging task, given the rich nonlinear dynamics the device is capable of showing and the strength of the nonlinear response. To accomplish these tasks, the oscillator is mounted on a shaking table and it is driven through several levels of excitations with both harmonic and random inputs.
e nonlinear system identification is performed by applying the nonlinear subspace identification (NSI) method [13][14][15]. Generally, NSI estimates a nonlinear state-space model starting from the input-output data, by separating the nonlinear part of the response from the socalled underlying linear system. is is not a straightforward operation for the case considered here though. e bistable nature of the device implies that the underlying linear system is not stable, as it exhibits a negative linear stiffness. erefore, NSI cannot be directly applied and some prior data manipulation is required. Nonetheless, the identification strategy proposed in this work is capable of estimating both the coefficients of the nonlinearities and the underlying linear systems associated with the stable equilibria starting just from one cross-well measurement, with satisfying accuracy. e nonexistence of a stable underlying linear system and the strong nonlinear behavior make the test set really interesting for the purposes of nonlinear system identification, confirming the effectiveness of the nonlinear subspace identification method. Eventually, the restoring force surface (RFS) method [16,17] is also applied as a comparison.

A Negative Stiffness Oscillator
e device here studied consists in a U-shaped steel frame connected through rods to a central moving mass. e frame has the purpose to keep the rods under compression during their movement so that a bistable mechanism is achieved. A schematic representation of the device can be seen in Figure 1. e lower surface of the frame is attached to a shaking table so that a displacement b(t) can be imposed to the structure. It is also assumed that the inertia of the moving parts can be concentrated into one central point with mass m, comprising the mass of the central bushing and the equivalent inertia of the rods. e vertical movement of this point is described by the coordinate y(t), while the rotation of the rods is called θ. e elasticity of the frame is also taken into account, and it acts like a compression force p(θ) to the rods. Since the flexural elasticity of the frame is much higher than the axial elasticity of the rods, the latter is considered infinitely rigid. A schematic representation of the functional model here described is reported in Figure 2, together with the free-body diagram of the mass m. e equilibrium equation along the vertical coordinate is m € y(t) + 2p(θ(t))sin(θ(t)) + mg � 0. (1) e system is a single-degree-of-freedom system; thus, just one variable is needed to describe the motion. In particular, since vertical displacements and accelerations are measured in the experimental setup, the variable z(t) is taken as an independent variable; therefore, both θ and p(θ) should be written as a function of z. e elasticity of the frame can be analyzed to obtain the force that the frame transmits to the rods. Just half of the frame is considered in the following, as depicted in Figure 3.
When the mass moves along the vertical axis, the frame bends and deforms until a new equilibrium position is reached, corresponding to the resting (undeformed) position of the frame. e analytic expression of p(θ) can be found by studying the flexibility of the half-frame, considered as two connected cantilever beams under bending stress. It is not the purpose of this work to analytically derive the expression of p(θ). Instead, a qualitative representation of its total vertical component p ζ (θ) � 2p(θ)sin(θ) is depicted in Figure 4 as a function of z(θ).
It can be seen that p ζ (z) has three roots and crosses the origin with a negative slope. Also, it can be approximated as a polynomial expansion of degree 3:  Shock and Vibration e equation of motion can eventually be written as Equation (4) has the form of a negative stiffness Duffing equation, whose coefficients k 3 , k 2 , and k 1 have to be estimated. In particular, the final objective of this work is to identify them directly from the experimental measurements, in order to assure a good correspondence between the numerical model and the real measured behavior. Also, the restoring force K(z) and the potential U(z) can be defined as A qualitative representation of the potential is shown in Figure 5, where its double-well nature can be clearly identified.
e potential is not symmetric because of the gravitational contribution and the asymmetry of the frame. Also, the three equilibrium positions are depicted, obtained by setting K(z * ) � 0. Two out of three positions represent a stable equilibrium, namely, z * − and z * + , while the central position z * 0 is an unstable equilibrium point.
When the oscillations are bounded around one equilibrium position z * ± , the system acts like a stable nonlinear oscillator, and the associated linear natural frequency ω ± can be computed by with U ″ (z * ± ) being the second derivative of U(z) computed in z * − or z * + . It is worth writing the equation of motion considering the oscillations around one of the possible equilibrium points. To ease the notation, the generic equilibrium position is called z * . A new variable x(t) can be defined as Substituting equation (8) in equation (4) yields with Finally, equation (9) becomes e advantage of using this formulation instead of equation (4) is that it allows the definition of a stable underlying linear system. is is a crucial requirement for the nonlinear subspace identification method, adopted in Section 4 to perform the nonlinear system identification of the structure under test.

Experimental Characterization
Two photos of the experimental setup corresponding to the two stable equilibrium positions are reported in Figure 6. e moving mass is instrumented with an accelerometer to measure its absolute acceleration € y(t) and a laser vibrometer to measure its absolute displacement y(t). e zero position of y(t) corresponds to the horizontal configuration of the rods (θ � 0). e acceleration of the base € b(t) is also recorded through a second accelerometer, and the displacement z(t) is computed as the difference between the laser measure y(t) and the displacement of the base b(t), obtained by integrating twice its measured acceleration.   Shock and Vibration e system is excited with both harmonic and random inputs to characterize its linear and nonlinear behavior.

Random Tests.
A first set of measurements is performed with a random input over the frequency range 7-50 Hz. e sampling frequency is f s � 512 Hz, and the duration of the acquisition is t � 300 s. Several forcing levels are applied, expressed as RMS values of the acceleration of the base € b, and the starting position is z * − . e results are depicted in Figure 7 in terms of time series and transmissibility T, defined as the ratio between the output acceleration € z and the input € b. It can be noted that the mean value of the displacements is not null, as in-well oscillations take place in the neighborhood of the equilibrium position, which in this case is approximately at −0.03 m from the horizontal configuration of the rods (θ � 0). Also, the asymmetry of the displacements increases as the forcing level, as the system tries to cross the negative stiffness region and to reach the positive equilibrium position.
is results in a clear change in the transmissibility, where a softening effect can be seen when increasing the excitation level, in accordance with the theoretical studies that show a similar behavior in the case of in-well motion [9].
If the energy given to the system is high enough, cross-well oscillations are retrieved and the moving mass oscillates in a wider region, including the two stable equilibrium positions z * − and z * + . is situation is represented in Figure 8, where the displacement clearly shows repeated crossings between negative and positive values. Also, the statistical distribution of z(t) depicted in Figure 8(b) highlights the asymmetry of the structure so that oscillations around the negative position are roughly the 66% of the total acquisition length. is result agrees with the shape of the potential retrieved from the model described in the previous section, which shows two wells with different heights. e final confirmation will be given in the following section where the experimental data are processed so as to identify the actual potential of the system.

Sine-Sweep and Harmonic Tests.
Harmonic excitations are a powerful tool to understand the behavior of nonlinear systems. is is particularly true for the case considered here, because of the possibility of chaotic motion. A logarithmic sine-sweep in the frequency range 7-21 Hz is considered first.
e sampling frequency is f s � 512 Hz, and the length of the up-down sweep is t � 240 s. Two forcing levels are taken into account, expressed as the amplitude b 0 of the base displacement, and the starting position is the negative equilibrium z * − . e up and down sweeps are shown in Figure 9. Both the excitation levels clearly produce a nonlinear behavior for the in-well motion, as classic jump phenomena can be observed in the up and down sweeps. e softening effect can be seen in Figure 9(a), where an increase in the forcing amplitude corresponds to an earlier occurrence of the jump-up. Also, two distinct jumps can be noticed, corresponding to the dominant frequency (10-11 Hz) and its second harmonic (20)(21).
If the input is high enough, cross-well motion is retrieved also in this case (Figure 9, orange line). It is interesting to look at the harmonic contributions in this case by computing the spectrogram of the relative displacement. e result is reported in Figure 10, where the first two minutes refer to the sweep up, while the second two minutes refer to the sweep down.
Both even and odd harmonics of the instantaneous frequency are present along the whole acquisition, confirming the asymmetrical behavior of the nonlinear system. Subharmonics are also visible in some regions, in particular around 2 minutes. Generally, they are symptomatic of the possibility of bifurcations and chaotic motion [18,19], and thus, a series of harmonic tests with constant frequency is performed to analyze these effects. e excitation frequency is ] � 9 Hz, and three different amplitudes b 0 are considered.
e results are presented in the phase diagrams in Figure 11. e different excitation amplitudes result in different kinds of motion of the moving mass, ranging from harmonic oscillations to cross-well motion. In particular, when the amplitude is b 0 � 2 mm (Figure 11(a)), the phase plane shows one closed orbit centered around the equilibrium position, i.e., one periodic solution [20]. When the  Shock and Vibration 5 amplitude increases to b 0 � 4.7 mm, some nested orbits can be seen in the phase diagram ( Figure 11(b)). Two paths in a closed loop are generally representative of the period doubling effect [7] so that a subharmonic with twice the period of the dominant shows up. Such a behavior is retrieved in the experimental case of Figure 11(b), and this is clearer when looking at the power spectral density Z of the displacement, represented in Figure 12. It can be seen that the system responds at both integer multiples of the dominant frequency ] (2], 3], . . .) and of its subharmonic (1/2)v ((3/2)v, (5/2)v, . . .). As for the highest excitation level in Figure 11(c), a crosswell motion occurs, and the solution is not periodic anymore.
is behavior lasts for the whole acquisition (10 minutes), and a portion of the time response can be seen in Figure 13.
is kind of persisting nonperiodic response to a periodic excitation is a symptom of chaotic behavior. It should be recalled that no definition of chaos is universally accepted, and this is particularly true when experimental data are considered. e reason is that uncertainties and noise in the data acquisition may interfere with the extreme sensitivity to the initial conditions that characterize chaotic systems. An useful way to check whether a system is behaving chaotically or not is to look at its largest Lyapunov exponent λ [20]. A positive sign of λ means chaotic motion, while a negative sign is a representative of a periodic orbit. Several methods exist to compute λ from experimental time series, and the one proposed in [21] is adopted here. e results are shown in Figure 14, where a positive λ is retrieved.
Eventually, the experimental Poincaré sections are computed for different phase synchronizations ϕ of the data with the forcing term [22]. e typical shape of a strange attractor is retrieved [8] and depicted in Figure 15(a) in a polar plot, while one of its sections is represented in Figure 15(b). e considered structure is then capable of exhibiting a variety of motions, as expected from the theory. Among them, the cross-well is surely the richest in terms of dynamics, as it covers all the possible positions of the moving mass. For this reason, cross-well measurements will be used in the following section to identify the nonlinear model parameters.

Nonlinear System Identification
A first guess of the nonlinear restoring force K and the potential U is estimated from the measured time series using the restoring force surface method (RFS) [16,17]. is method is fairly simple and allows to visualize the nonlinearity, providing a very useful insight to the user. On the contrary, its capabilities are very limited when trying to identify a nonlinear model structure, as it is essentially based on raw data processing and simple fitting. For this reason, the final model is identified using the nonlinear subspace identification method (NSI) [13][14][15]. NSI gives as outcomes the fully nonlinear state-space representation of the system, together with the FRF of the so-called underlying linear system and the estimation of the coefficients defining the nonlinearities. It requires the input-output data and the knowledge of the nonlinear basis functions, whose selection can take advantage from the preliminary RFS analysis.

First Step: RFS Method.
e equation of motion describing the dynamical system considered here can be written in the following form: where f(t) is the forcing term and R(z, _ z) is the nonlinear restoring surface, generally a function of the displacement z and the velocity _ z. If the inertial term is shifted to the righthand side of the equation, the restoring surface can then easily be computed and its features extracted. In particular, if small velocities are taken into account, such that | _ z| < ε s , the obtained slice of the restoring surface approximates the restoring force K(z). On the contrary, when small displacements around the equilibrium positions are considered, such that |z − z * | < ε d , an approximation of the damping force, called D( _ z), is retrieved. For the case considered here, the cross-well measurements are used to build R(z, _ z). e velocity _ z is obtained by integrating the accelerations € y and € b and subtracting them.
e obtained restoring surface is reported in Figure 16, together its sections K(z) and D( _ z), computed setting ε s � ε d � 0.1%. Shock and Vibration e restoring force is eventually fitted to a polynomial expansion of degree 3 as in equation (5), obtaining where k 0 is the static component due to gravitational effects. Once the expression of K(z) is obtained, the potential U(z) can be computed from equation (6). Both K(z) and U(z) are shown in Figure 17. As for the theoretical model presented in Section 2, the potential shows two wells with different heights and three equilibrium positions are identified, two stable and one unstable. e final values of the coefficients defining the nonlinearity are estimated using NSI, and a comparison with the ones obtained using RFS is also reported in Section 4.2.
As for the damping plot in Figure 16, it is rather difficult to estimate a proper damping model due to the very high dispersion of D( _ z). Nevertheless, a possible and realistic guess is that some friction is present between the bushing of the moving mass and the vertical steel guide. For this reason, a nonlinear damping function is also added to the set of nonlinear basis functions given to NSI in the following.

Second
Step: NSI Method. e nonlinear subspace identification method is based on a nonlinear state-space representation of the system, obtained by considering the nonlinear terms as feedbacks to the underlying linear system and processing the measured input-output data using the subspace formulation. e existence of the underlying linear system is an essential requirement for the method to work. In the present study, this requirement is not fulfilled when the system goes through cross-well oscillations because of the bistable nature of the device. A possible overcome to this  Shock and Vibration issue has been already shown in Section 2, and it consists of a simple shift of the reference axis with respect to one stable equilibrium position z * , called "reference value" hereafter. It should be noted that the two stable equilibrium positions of the system correspond to two underlying linear systems mutually exclusive, each one existing only when the moving mass oscillates around that equilibrium position. us, the two underlying linear systems can be identified with NSI from a single cross-well measurement by selecting the corresponding equilibrium position as the reference value. Whatever reference position is used, a new displacement variable x(t) � z(t) − z * can be defined, as in equation (11). e coefficients of the two nonlinear functions become k 3 and k 2 , and nonlinear damping is also considered in the form _ x| _ x|. Equation (11) can then be rewritten according to the feedback formulation: where a linear viscous damping contribution is also considered and R nl (x, _ x) is the assumed nonlinear restoring surface, equal to Referring to Section 3.1, the random test with excitation amplitude 38 m/s 2 RMS (Figure 8) is considered for the nonlinear system identification with NSI. In particular, the last 60 seconds are used as a validation set for the evaluation of the residuals with the measured output, while the identification is performed on the remaining part of the acquisition. e inertance of the system G is depicted in Figure 18, together with the coherence estimation. As expected, the FRF is very distorted due to the strong nonlinear contributions, and the coherence drops to very low values in the frequency region up to the resonance peak.
NSI is performed considering first the negative equilibrium position as a reference value, meaning z * � z * − . Stability is checked varying the model order for frequencies, damping ratios, MACs, and modal masses [15], and the stabilization diagram is reported in Figure 19. Since the modal parameters show stability from a model of order 2, this is the selected value.
Once the model order has been selected, the fully nonlinear state-space model is retrieved, together with the underlying linear system and the coefficients of the nonlinear basis functions. is is repeated also when the reference value is equal to the positive equilibrium position, meaning z * � z * + . e stabilization diagram for this case is not reported since it is similar to the one in Figure 19, and a model order equal to 2 is accomplished also in this case. e identified modal parameters of the two underlying linear systems related to the two reference values are reported in Table 1 in terms of natural frequencies and damping ratios.
e FRFs of the two underlying linear systems are depicted in Figure 20 together with the measured (nonlinear) one, already shown in Figure 18. e coefficients of the nonlinear basis functions can also be computed from the nonlinear state-space model, recalling that they are treated by the method as frequency-dependent and complex-valued quantities [13]. Assuming that the true coefficients are real numbers, the imaginary part of the identified counterparts should then be zero and the real part should not depend on the frequency. Since this happens only in complete absence of noise and nonlinear modeling errors, the ratio between real and imaginary parts can be taken as an indicator of the goodness of the nonlinear basis functions choice.
Running NSI for the two different reference values means that not only two independent underlying linear systems are retrieved but also two nonlinear state-space models. is result in a double estimation of the coefficients of each nonlinearity. In particular, the estimation of the    coefficients k 3 and c nl should be the same when NSI is applied to the two reference positions, as they are both invariant in the equation of motion (equations (14) and (15)). On the contrary, k 2 depends on the choice of z * . e identified coefficients are depicted in Figure 21 in their real parts for the two reference positions. A summary of the identified coefficients is also reported in Table 2 for the two reference values. e final value of each coefficient is computed as the spectral mean of its real part over the selected frequency range, and the average ratio between real and imaginary parts E[R/I] is also shown.
It is worth highlighting that the imaginary part is always much lower than the absolute value of the real part in the selected frequency range, which assesses the goodness of the identification. Also, the real parts of k 3 and c nl computed with respect to the two reference values are very close to each other, as expected. Since the system has a "favorite" equilibrium position, which is the negative one, the number of samples associated with oscillations around this position are higher than the other. For this reason, the estimation of the "negative" underlying linear system is considered more reliable and so is the identification of the nonlinearity. e  chosen final values of the coefficients are then the ones related to the negative reference position. ese coefficients can then be compared with the ones estimated by the RFS method, by looking at the identified restoring force. e comparison is reported in Figure 22, where the damping force estimation is also shown. e agreement between RFS and NSI estimation of K(z) is very high so that the two curves in Figure 22(a) appear to be almost overlapped, and the percentage deviation is around 1%. As for the damping force D( _ z), it is difficult to evaluate the goodness of the estimation due to the high dispersion of the points. However, at least for the linear   part (for small | _ z|), it is possible to see that the slope of the NSI estimation and the data points agree.
Finally, the identified nonlinear state-space model can be used for reproducing the output time series, given the input force. Calling the simulated output x NSI , a comparison with the measured one can be carried out to validate once more the identification. e validation set corresponds to the last 60 seconds of the acquisition, and the results are reported in Figure 23 as simulated time history against its residual with the measured one in time and frequency domains. e identified state-space model is capable of catching the cross-well motion with a very good accuracy, providing a percentage RMS deviation from the measurement of approximately 8%.

Conclusions
A negative stiffness oscillator has been characterized experimentally to exploit its dynamical behavior. A variety of different kinds of motions can be obtained because of the bistable nature of the device, from in-well to cross-well oscillations, including chaotic motion. ese dynamical behaviors have been confirmed by the experimental observations, retrieving linear, nonlinear, and chaotic motions. Eventually, the parameters defining the nonlinearity have been recognized via nonlinear system identification, adopting two different methods. A first guess has been obtained using the restoring force surface method, whose implementation allows to easily visualize the nonlinear behavior and the asymmetric double-well potential of the system. e final identification has been performed using the nonlinear subspace identification (NSI) method with a cross-well random measurement. NSI proved to be a robust technique, capable of handling a very strong nonlinear behavior. e nonlinear coefficients have been eventually estimated, together with the linear modal parameters associated with small oscillations around the two equilibrium positions.

Data Availability
e data used to support the findings of this study are available from the corresponding author upon request.