Experimental Characterization of Friction in a Negative Stiffness Nonlinear Oscillator

Nonlinear dissipative phenomena are common features of many dynamical systems and engineering applications, and their experimental characterization has always been a challenge among the research community. Within the wide range of nonlinear damping mechanisms, friction is surely one of the most common, and with a high impact on the dynamical behavior of structures. In this paper, the nonlinear identification of friction in a negative stiffness oscillator is pursued. The structure exhibits a strong nonlinear behavior, mainly due to its polynomial elastic restoring force with a negative stiffness region. This leads to an asymmetric double-well potential with two stable equilibrium positions, and the possibility of switching between them in a chaotic way. Friction plays a crucial role in this context, as it derives from the continuous sliding between the central guide and the moving mass. The system is driven through harmonic tests with several input amplitudes, in order to estimate the variations in the energy dissipated per cycle. The identification of the frictional behavior is then pursed by minimizing the errors between the experimental measurements and the model predictions, using the harmonic balance method in conjunction with a continuation technique on the forcing amplitudes.


Introduction
Friction is a complex nonlinear phenomenon highly present in mechanical systems. It results from the interaction between adjacent surfaces and is dependent on their topography and materials, and on the presence of lubrication and relative motion [1]. Studies on friction have been carried out for more than 500 years [2], and a more detailed explanation about modeling techniques can be found in several works in the literature [3][4][5]. A major classification of friction is usually done by separating the "pre-sliding" from the "sliding" (or "gross-sliding") behavior. The first can be associated with the elastic and plastic deformations at the asperity levels of the adjacent surfaces, and it is known to be mostly dependent on their relative displacement. The friction force in the pre-sliding regime is a hysteretic function of the position with a non-local memory [1,4]. Instead, the sliding regime is largely due to the shearing resistance of the asperities, and mainly dependent on the velocity. The pioneering friction model of Coulomb [6] is indeed the most known example of sliding friction.
The experimental identification of friction has always been a difficult task, especially when reallife mechanical structures are involved. This is indeed due to the intrinsic nonlinear nature of the friction phenomenon, which may lead to limit cycles, tracking errors, stick-slip motion, hysteresis and other typical features of nonlinear systems such as natural frequency shifts and jumps. Many nonlinear identification techniques have been developed by the research community in the last decades, and two extensive literature reviews can be found in [7,8]. Nevertheless, different methods are suitable for the different classes of nonlinearities and kind of excitations, and few of them have been tested with frictional and/or hysteresis nonlinearities. In [1], several identification techniques of pre-sliding and sliding friction were compared, ranging from white-box to black-box models. The identifications were performed directly on the friction force, which was measured in a controlled laboratory setup with a tribometer. In [9], the same experimental setup was adopted and the identification was performed considering the LuGre [10] and the Maxwell-slip [11] friction models in conjunction with nonlinear regression techniques. The LuGre model was again adopted in [12] to identify the frictional behavior in a torque motor control system, using an evolutionary algorithm to minimize the residuals between predictions and measurements with constant velocity experiments. A different approach was adopted in [13], where a general hysteresis model was identified via the nonlinear state-space modeling of a single-degree-of-freedom system under multisine excitation. Ground vibrations sine sweep tests were conducted in [14], to estimate the shape of the friction force at the joints of the wing-payload substructure of an F-16 aircraft using the restoring force surface (RFS) [15] and wavelet transform methods. Indeed, friction phenomena associated with joints are known to be a typical source of nonlinearity in mechanical systems [16], and this represents a critical issue in several engineering applications, such as bladed disks [17]. In such cases, friction is usually modeled with node-to-node contact elements, where Coulomb's friction law is used to model the local slip conditions of the joint interfaces [18]. On the experimental side, the nonlinear dynamical behavior of contact interfaces has been widely observed, for instance in [19][20][21].
In this paper, the data-driven nonlinear modeling of the frictional behavior of a negative stiffness oscillator is performed. The 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 reduce the wave propagation, known to be a critical phenomenon in such structures [22,23]. It is characterized by a strong nonlinear behavior mainly due to its double-well characteristics [24], and it exhibits two stable equilibrium positions plus an unstable one. The oscillations can either be bounded around one stable point ("in-well") or include all the three positions ("cross-well"). In both cases, periodic oscillations can evolve to steady in-well or cross-well chaotic motions under external periodic excitations [25,26]. The bi-stable nature of the device is translated to a polynomial-kind restoring force with a negative slope around the origin. This has been experimentally identified in a previous work [27] using the nonlinear subspace identification method (NSI) [28] with cross-well oscillations under random excitation. However, the method struggled when trying to infer a nonlinear dissipative model.
The identification of the friction force of the device turns out to be a particularly challenging task, as it is not the only and dominant source of nonlinearity. The device is here driven through several harmonic excitations with increasing amplitudes, and the friction force is identified by minimizing the residuals between the model predictions and experimental measurements of the energy dissipated per cycle over the considered input amplitude range. The model response is computed using the harmonic balance method (HBM) [29], therefore searching for periodic responses approximated by Fourier series. HBM is usually adopted in conjunction with a continuation technique [30,31] to study the evolution of the periodic solutions with respect to the excitation frequency and to build the so-called "nonlinear frequency response curves". The same idea is used in this paper, but a novel procedure is implemented adopting the continuation technique with respect to the excitation amplitude. This allows to estimate the "nonlinear amplitude response curves", and therefore the dissipated energy per cycle. Given the interest in obtaining a physically based model, a white-box approach is pursued in this work. A proper frictional model should then be considered, which is first inferred using the RFS method together with physical insights about the device functionality. A genetic algorithm [32] is eventually used to optimize the model parameters based on the experimental observations.
Results show a good confidence in the identified parameters and provide a reliable model able to catch the strong nonlinear dynamics of the structure under test. This confirms the effectiveness of the presented methodology, which can be applied to a wide class of nonlinear systems.

A Negative Stiffness Oscillator
The device consists of a U-shaped frame connected through rods to a central moving mass. The frame keeps the rods under compression during their movement, achieving a bi-stable mechanism. A schematic representation of the device is depicted in Figure 1. The lower surface of the frame is attached to a shaking table, so as to impose a displacement ( ) to the structure. An accelerometer is located on the shaking table to measure the acceleration of the base ̈( ). The displacement ( ) is then obtained by integrating twice its measured acceleration. It is also assumed that the inertia of the moving parts can be concentrated into one central point with a mass of = 0.26 , comprising the mass of the central bushing and the equivalent inertia of the rods. The vertical movement of this point is described by the coordinate ( ) and it is measured by a laser vibrometer. The zero position of ( ) corresponds to the horizontal configuration of the rods. The reader can refer to [27] for more details about the device.
It is assumed that the system can be modeled as a single-degree-of-freedom system in the variable ( ) = ( ) − ( ), leading to an equation of motion of the kind ̈+ ℛ( ,) = −̈= ( ), (1) where ℛ( ,) is the restoring surface, the sum of the (elastic) restoring force ( ) and the damping force ( ,). The latter is supposed to be a function of both displacement and velocity ̇. The restoring force of the considered system can be expressed as a polynomial expansion of degree three: Therefore, the undamped equation of motion takes the form of a double-well asymmetric Duffing oscillator [26], and its potential can be defined as A qualitative representation of the potential is shown in Figure 2, where its double-well nature can be clearly noticed. The potential is not symmetric because of the gravitational contribution and the asymmetry of the frame. Further, the three equilibrium positions are depicted, obtained by setting ( * ) = 0. Two out of three positions represent a stable equilibrium, namely − * and + * , while the central position 0 * is an unstable equilibrium point. The oscillations of the moving point are said to be "in-well" when the motion is bounded around one of the two stable equilibrium positions ± * . The associated linear natural frequency ± can be computed by ′′� ± * � being the second derivative of ( ) computed in − * or + * . Two photos of the device in correspondence with the two stable equilibrium positions are depicted in Figure 3.

Experimental Nonlinear Characterization and Equation of Motion
Both random and harmonic tests have been performed to identify the features of the device, and all the tests have been conducted with a sampling frequency of 512 . In particular, random tests have been used in a prior work to identify the restoring force of the device, whose results are briefly summarized in the following. Harmonic tests are instead used in this work to identify the damping force. The excitation force is expressed in this case in terms of excitation amplitude 0 and frequency .

Identification of the Restoring Force
The values of the coefficients of the restoring force have been experimentally identified in [27] using the nonlinear subspace identification (NSI) algorithm with cross-well measurements under random excitation. A portion of the time history of the measured displacement ( ) is reported in Figure 4, while the identified values are reported in Table 1.
According to Equation (4), the linear natural frequencies of the small oscillations around the two stable equilibrium positions are − = 11.6 ⋅ 2 / and + = 8.6 ⋅ 2 / . The RFS method can be applied to the time history of Figure 4 to visualize the experimental restoring and damping forces, and the results are depicted in Figure 5a. In particular, if small velocities are taken into account, such that || < , the obtained slice of the restoring surface approximates the restoring force . Figure 5b shows the experimental restoring force compared to the identified one of Table 1, which nicely fits the experimental observations. Instead, when small displacements around the equilibrium positions are considered, such that | − * | < , an approximation of the damping force should be retrieved. This is reported in Figure 5c, but the dispersion of the points appears to be too high to correctly deduce any damping model.

Characterization of the Friction Force
Given the difficulty of extrapolating a damping model from random measurements, two series of stepped sine tests have been performed with increasing forcing input amplitudes at a fixed excitation frequency of = 9 , starting from the two equilibrium positions. An extract of the recorded displacement is depicted in Figure 6.  | −| < . Results are represented in Figure 8, where different shapes of the damping force can be observed to be varying in position . The damping force generally resembles the sliding friction behavior, but it appears to vary with the position of the moving mass. In particular, the slopes of the negative and positive velocity regions change due to the polynomial restoring force. The amplitude of the damping force changes as well, and its variations can be described by a non-constant normal force along . Indeed, the normal force comes from the interaction between the moving mass and the rods, which varies during the movement of the mass. Given the above considerations, the damping force is assumed as follows: Equation (5) represents a continuous friction model, as in [33], and includes the following: • A dynamical damping friction term ( ), with a Coulomb-like function having a transition velocity ; • A static friction term ( ) to account for the stiction force and the Stribeck effect; • A linear viscous damping term ̇, to account for possible viscous forces generated by the lubricated slider.
The dynamic friction term ( ) is assumed to be non-negative and described by a quadratic function, while the static term ( ) is assumed to be proportional to the dynamic term, yielding The complete equation of motion of the system is therefore The unknowns to be identified are the parameters defining the damping force, i.e., , , , , and .

Optimization-Based Identification Strategy
The identification pursued in this work aims to estimate the damping force via an optimizationbased strategy. In this context, the definition of the cost function to be minimized is a crucial step, as it has a huge influence on the feasibility of the optimization. This is usually chosen as the residual between a measured feature and the corresponding model prediction. A valid feature in this context is the energy dissipated per cycle (EDC), given by the following considerations: • The external forcing term is periodic, and the response is expected to be periodic as well (excluding the chance of chaotic motion); • The evolution of the dissipated energy with the input amplitude is directly correlated with the dissipation functional form, i.e., the damping force ; • It depends both on the response amplitude and the phase with the forcing term.
The ℎ measured EDC corresponding to an input amplitude 0, of the stepped sine test is called , while the corresponding model prediction is called . The cost function is therefore defined as the sum of the absolute differences: Model predictions depend on the set of parameters to be optimized, which can be recast into a vector as follows: The final value � is the minimizer of the cost function ( ): � = arg min ( ). (10) Note that the vector of parameters takes into account the values of the dynamic friction forces in − * , + * , 0 rather than , , of Equation (6) to have a better physical interpretation of the outcome. The coefficients , , can eventually be found by fitting the quadratic function ( ) to the considered three points.
A genetic algorithm is adopted in this work to find the best set of parameters � . Genetic algorithms belong to the class of evolutionary global optimizers, and they are commonly used to generate high-quality solutions to optimization problems using biologically inspired mechanisms, such as reproduction, mutation and selection. Candidate solutions act like individuals in a population, which evolves through successive generations. A portion of the existing population is selected at each generation to breed a new offspring, and the selection is made upon the corresponding values of the cost function. Modifications can be introduced to better explore the range of possible solutions and avoid local minima. For instance, a mutation rate can be defined to introduce random changes to the existing solutions. The reader can refer to [32] for an exhaustive description of evolutionary (and genetic) algorithms.

Experimental Estimation of the Energy Dissipated per Cycle (EDC)
The EDCs are directly estimated from the experimental measurements by numerically computing the integral of the response ( ) with respect to the harmonic input ( ) for each cycle and averaging over the number of cycles. Transients are removed, so as to take into account only the steady-state responses. Since the measurements are performed starting from the two equilibrium positions, two sets of dissipated energies are collected, called − and + and refer to the oscillations around − * and + * , respectively. The complete vector of measured EDCs is therefore = vec{ − , + }. Each final value is obtained by averaging over the cycles, and the associated standard deviation is used as an index of dispersion.
A representative response-input plot is reported in Figure 9 with an excitation amplitude equal to 1.3 , while Figure 10 depicts the evolution of with respect to the input amplitudes vector 0 . The range of considered input amplitudes goes from 0.2 to 2.7 . Period doublings start to occur in the experimental measurements beyond this value, which are not the object of this study. The error bars in Figure 10 represent the quantity ± 3 .

Numerical Computation via Harmonic Balance Method with Amplitude Continuation
Considering a harmonic input of the kind ( ) = 0 sin( ), the energy dissipated over one cycle can be written as Since the system described by Equation (7) is nonlinear, the response ( ) to a periodic excitation can be written in terms of a Fourier series with H harmonics: �. (12) It can be easily proved that the integral of Equation (11) is zero for each integer harmonic contribution of the response different from 1. Therefore, the EDC reduces to (13) where 1 is the phase of the first harmonic of the response. Note that Equation (13) is valid only if sub-harmonics are not taken into account.
The Fourier coefficients of the response ( ) can be rapidly computed via the harmonic balance method for each forcing input amplitude 0, and for a fixed . In particular, the focus is to track the evolution of varying 0 , and to compare the model predictions with the experimental observations of Figure 10. An arc-length continuation procedure (like the one proposed in [30]) is adopted to study the evolution of the dissipated energy. The forcing input amplitude 0 is therefore chosen as the "continuation parameter" and the solution for each 0 is sought along the arclength of a solution branch. This allows to account for possible singularities of the system being Jacobian, which lead to response bifurcations [29]. An illustrative example is exposed in the following.

Illustrative Example
A linear oscillator with continuous velocity-based friction is considered in this example. The system parameters are summarized in Table 2 and the equation of motion reads The excitation frequency is set to 80% of the linear natural frequency of the system � / , and the amplitude increases in the range 0 − 2 . Seven harmonics are included in the response. The nonlinear amplitude response curve obtained using HBM with the amplitude continuation is depicted in Figure 11, together with the EDCs and the evolution of the harmonic coefficients. An unstable branch can be noted in the region between 1.3 and 1.4 . Multiple solutions exist in this region, meaning that the system response would avoid the unstable branch and suddenly reach a new stable solution. This phenomenon is very common in nonlinear systems when sweeping around their resonance frequencies and it is usually called "jump" [29]. Here, the nature of the phenomenon is different, but the behavior is very similar. Figure 12 depicts the response of the system around the unstable region computed by numerically integrating the equation of motion with the Newmark algorithm [34] and under the amplitude sweep up and down excitations. The nonlinear amplitude response curve computed with the HBM is also overlapped. The jump phenomenon can be clearly observed. A slight difference can be noted between the HBM and numerical simulations around the jump region. This is probably due to the different kind of excitation in the two approaches: numerical simulations are conducted with an amplitude sweep excitation, which contains transients. The HBM represents instead the steady-state solutions for each input amplitude value.

Results and Discussion
The unknown parameters of Equation (7) are identified by optimizing the residuals between the measured EDCs of Figure 10 and the predicted ones according to Equation (13). Both sets of measurements around the positive and negative equilibrium positions are considered.
The genetic algorithm is applied considering a population of 100 individuals, 20 generations and 20% mutation rate. The latter in particular is to allow the algorithm to explore a wider range of possible solutions, by randomly modifying a portion of the existing population. The cost function with respect to the single parameters of the optimization is depicted in Figure 13. It can be noted that the minimum value of the cost function is reached sharply when considering the parameters ( − * ), (0), ( + * ) and , while results are more dispersive for the others. It should be highlighted though that these four parameters have the biggest influence in the EDC of the considered measurements. The transition velocity and the static proportionality coefficient have a major influence in the low-amplitude region, and tend to compensate each other. Some more details are presented in the sensitivity analysis in the next section.
The value of the cost function across the generations of the genetic optimization is reported in Figure 14. The minimum value of the cost function is ( � ) = 0.012, and the corresponding set of parameters � is listed in Table 3.  Eventually, the comparison between the measured EDCs and the final model predictions is reported in Figure 15. The average percentage deviation between the measured and predicted EDCs is 6% for the positive oscillations and 12% for the negative ones. The model nicely captures the nonlinear dynamics of the system, especially at higher input amplitudes. Errors are instead more consistent in the low-amplitude region. Nevertheless, displacements in this region are so small that other phenomena might be involved, such as displacement-dependent and hysteretic frictional behavior.

Sensitivity Analysis
A sensitivity analysis is carried out on the parameters of the identification to check their influence on the predicted EDCs. The parameters are varied in the measure of 20%, 50%, 70%, 120% and 150% of their identified values in Table 3. In particular, ( ) is varied with a coefficient , so as to consider the quantity ⋅ ( ). Figure 16 shows the results of the sensitivity analysis on the EDCs associated with the oscillations around the positive equilibrium position. dispersion in Figure 13. Furthermore, low-amplitude oscillations might be affected by other phenomena, such as hysteresis loops, which are not considered in this study.

Nonlinear Amplitude Response Curves
Finally, the identified model is used to build the nonlinear amplitude response curve of the structure, with an excitation frequency = 9 and an input amplitude range of 0.1 − 3.7 . As already pointed out in Section 4.1, period doublings start to occur in the experimental measurements after 2.7 . For this reason, sub-harmonics up to 1/4 are also included in the following test to check whether this happens also in the model. Figure 17 depicts the nonlinear amplitude response curve starting from the positive equilibrium position, as well as the normalized harmonic coefficients. It can be noted that sub-harmonics start to show up around the 3 input amplitude, resulting in a period-doubling cascade. Representations like the one in Figure 17a are usually called "bifurcation maps". Each point of a bifurcation map represents the amplitude(s) of the steady-state solution for a specific value of the excitation amplitude. However, the map in Figure 17a is built using the harmonic balance method, which allows only a certain number of sub/super-harmonics to be considered in the response. In other words, it is not possible to predict chaotic behavior with this tool, as a periodic response is needed. A final validation of the goodness of the model is depicted in Figure 18, showing the experimental bifurcation map compared with the predicted one. The deviations between the predictions and observations resemble the ones obtained with the EDCs in the amplitude range up to 2.7 N. Higher discrepancies show up after period doublings start to appear. This is somehow expected, since the model has been trained on lower amplitudes. Further, period doublings are difficult to handle experimentally, as the system starts to show a considerable dependence on small external perturbations, such as noise. This makes the amplitudes of the periodic solutions to be quite variable, and the forcing input difficult to control. The encouraging results prove the goodness of the presented methodology and provide a reliable model for a system with rich and strong nonlinear dynamics.

Conclusions
In this paper, the nonlinear identification of a negative stiffness oscillator has been presented. In particular, the damping force of the device has been identified upon a physically based model, governed by a frictional behavior and dependent on the velocity and position. The identification has been carried out by minimizing the residuals between the model predictions and experimental observations. In particular, the energy dissipated per cycle has been adopted as the optimization criterion. On the experimental side, the system has been driven through harmonic tests with several input amplitudes, in order to estimate the variations in the energy dissipated per cycle. On the computational side, the harmonic balance method has been implemented with a continuation technique on the forcing input amplitudes. Results show a good confidence in the identified parameters, with low residuals in the analyzed input amplitude range. This confirm the effectiveness of the presented methodology, which can be applied to a wide class of nonlinear systems.