## Abstract

Strongly nonlinear structures have attracted a great deal of attention in energy harvesting and vibration isolation recently. However, it is challenging to accurately characterize the nonlinear restoring force using analytical modeling or cyclic loading tests in many realistic conditions due to the uncertainty of installation parameters or other constraints, including space size and dynamic disturbance. Therefore, a displacement-measurement restoring force surface identification approach is presented for obtaining the nonlinear restoring force. Widely known quasi-zero stiffness, bistable, and tristable structures are designed in a cantilever-beam system with coupled rotatable magnets to illustrate the strongly nonlinear properties in the application of energy harvesting and vibration isolation. Based on the derived physical model of the designed strongly nonlinear structures, the displacement-measurement restoring force surface identification with a least-squares parameter fitting is proposed to obtain the parameters of the nonlinear restoring force. The comparison between the acceleration integration and displacement differentiation methods for describing the restoring force surface of strongly nonlinear structures is discussed. Besides, the influence of the noise level on identification accuracy is investigated. In experimental conditions, quasi-zero stiffness, bistable, and tristable nonlinear structures with various geometrical parameters are utilized to analyze the identified nonlinear restoring force curve and measured force–displacement trajectory. Finally, experimental results verify the effectiveness of the displacement–measurement restoring force surface method to obtain the nonlinear restoring force.

## 1 Introduction

In recent years, many strongly nonlinear structures have been used in passive vibration isolation and energy harvesting. The basic principle of these nonlinear structures is to form an ideal nonlinear restoring force by coupling negative stiffness forces based on linear structures, such as an oblique spring/bar [1,2], X-shaped structures [3,4], cam-roller-spring mechanisms [5,6], a buckled beam [7,8], and magnetic coupling structures [9,10]. These nonlinear elements change the system's dynamic response and make it useful in many engineering applications, but they also increase the nonlinear restoring force calculation and identification complexity.

In the last decade, there were analytical calculations and cyclic loading testing methods to obtain the nonlinear restoring force. Carrella et al. [1] calculated the force–displacement function of the oblique spring-type high-static-low-dynamic-stiffness isolator and then simplified it with a cubic polynomial function. In complex spatial geometric structures such as diamond, X-shaped, and bionic structures invented by Jing and coauthors [3,4,11], the nonlinear restoring force of the dynamic equation was calculated using Lagrange's equations. Zhou et al. [5] utilized cam-roller-spring mechanisms to construct a quasi-zero stiffness force and replaced the exact expression of the nonlinear restoring force with an approximate piecewise form. These nonlinear restoring force calculation and simplification methods for vibration-isolation structures were always based on the physical model using ideal assumptions. For bistable or tristable structures with magnetic force coupling being utilized in energy harvesting, the magnetic dipoles theory, equivalent magnetic current, and magnetic charge methods were used to calculate the magnetic force [12–14]. However, these analytical approaches are limited by the magnet's dimensions and intervals, the accurate calculation of magnetic-field distribution along conductors. Moreover, Zhou and coauthors [9,10] and Leadenham and Erturk [15] measured the nonlinear restoring force of multistable beams and an M-shaped piezoelectric energy harvester and then fitted them with an empirical polynomial function for dynamic analysis. However, the estimated restoring force function cannot be adapted to parameter changes or if force-measuring instruments cannot be installed and arranged.

In addition to analytical calculations and direct measurement methods, identification strategies for the nonlinear restoring force, such as the nonlinear subspace method and restoring force surface method, had also received considerable interest in recent years [16]. For the nonlinear subspace-based method, a nonlinear restoring force was applied to an underlying linear system as a feedback force and then identifies the coefficients using a frequency-response function. Marchesiello and Garibaldi [17] numerically investigated this method to identify cubic and clearance-type nonlinearities. Different nonlinear frequency-domain subspace methods had been widely used in aerospace structures, such as impacts on mechanical stops [18], bolted connected joints on wingtips [19], and solar array structures [20]. However, the expression of the nonlinear feedback force requires a priori knowledge to ensure the accuracy of identification. The restoring force surface method, also called force-state mapping, was initially proposed by Masri and Caughey [21] and Crawley and Aubert [22] and had also been widely employed in nonlinear force identification. Al-Hadid and Wright [23] utilized the acceleration-measurement force-state mapping method to identify weakly nonlinear T-shaped beams and indicated that this method was more sensitive in the case of low damping. Worden [24,25] compared the integration and differential of measured time data to construct the restoring force surface and illustrated the influence of the selection of excitation signal on identifying the restoring force of linear and hardening spring oscillators. Kerschen and coauthors [26,27] employed the restoring force surface method to identify nonlinear piecewise beam and the helical wire rope isolator based on acceleration measurements. Allen et al. [28] proposed a piecewise-linear restoring force surface to identify the nonlinear electrostatic force of micro-cantilever beams based on velocity measurements. The wavelet transforms and restoring force surface method were combined by Noël et al. [29,30] to investigate the nonlinear behaviors of a real-life aerospace structure with multiple mechanical stops as well as F16 aircraft wing-to-payload mounting interfaces. Fuellekrug and Goege [31] studied the nonlinear restoring force identification of a large transport aircraft under two coupled modes. Moore et al. [32] used acceleration measurements to identify the nonlinear restoring force in a grounded nonlinear energy-sink system. The third-order Butterworth high-pass filter was introduced to reduce low-frequency components in the integration procedure and the trend term. It is found from the above work that the restoring force surface method and its improved version can identify the nonlinear restoring force without an exact analytical expression in realistic structures. Furthermore, the acceleration-measurement restoring force surface method has the advantages of low cost and easy implementation. However, when considering the identification of strongly nonlinear structures, especially those working in low-frequency vibration environments or with multistable properties, the previous acceleration and velocity measurement restoring force surface method will encounter the problem of low-frequency components in the integration procedure and complex trend items. Additionally, for flexible structures such as cantilever-beam systems, the influence of the size and mass of the additional acceleration sensors attached to the system must be considered carefully.

Therefore, the displacement-measurement restoring force surface approach is adopted in this paper to identify strongly nonlinear structures, and a comparison with the acceleration-measurement restoring force surface method is discussed. To design a test rig for accurate verification, quasi-zero stiffness, bistable, and tristable structures are implemented on a geometrically nonlinear cantilever beam with magnetic coupling. Based on the derived nonlinear model of the designed strongly nonlinear cantilever beam, the displacement-measurement restoring force surface identification method is adopted to obtain the nonlinear stiffness and damping force. The influence of noise intensity on identification accuracy is also investigated. Moreover, experimental measurements of quasi-zero stiffness, bistable, and tristable nonlinear structures with various geometrical parameters are performed to analyze the identified nonlinear restoring force and measured force–displacement trajectory. Experimental results verify the effectiveness of the proposed method.

The paper is organized as follows. In Sec. 2, the modeling theory of geometrically nonlinear cantilever beam is derived. The displacement-measurement restoring force surface identification process and comparison with acceleration-measurement restoring force surface method are explained in detail in Sec. 3. Section 4 is dedicated to the numerical investigation and comparison of the proposed method, and further experimental verification and evaluation can be found in Sec. 5. Finally, the essential conclusions of the present study are summarized in Sec. 6.

## 2 Modeling of Strongly Nonlinear Structures

A horizontal cantilever beam with a rotatable magnetic coupling is designed as shown in Fig. 1(a) to understand the dynamics of strongly nonlinear structures. A beam of length *L* carried two concentrated-tip magnets C with a mass of *M _{m}*. Magnets C are subjected to two rotatable magnets A and B, arranged opposite of magnets C. By changing the magnets’ interval or rotatable angles, the cantilever beam is made to exhibit various nonlinear characteristics [9,10]. In this paper, only the quasi-zero stiffness, bistable, and tristable beams are considered. Figure 1(b) shows a bistable stiffness beam under large deformation. The horizontal and vertical elastic displacements at the tip mass are

*x*and

*y*, respectively, and

*s*represents the distance along the beam's neutral axis. When considering an arbitrary point on the beam

*p*at a distance

*s*from the base, the laser displacement sensor can measure the response

*x*of the point.

_{p}*ρ*, cross-sectional area

*A,*Young's modulus

*E,*and second moment of inertia

*I.*Neglecting the effect of the rotary inertia of the beam and magnetic mass [33], the kinetic energy of the geometrically nonlinear beam-mass system is

*q*denotes the order of the polynomial function and

*k*is the coefficient. The relationship between the curvature and slope of the beam is

*φ*may be written in terms of the beam's elastic displacement as

_{p}*p*can be represented as a function of the mass tip displacement using the modal shape function,

*ψ*(

*s*), given as

*s*= 0, for example:

*x*, is obtained as

*N*

_{1}to

*N*

_{7}are given by

*Q*(

*t*) represents a generalized force, and a linear viscous damping force is considered in any part of the beam

It can be found from Eq. (16) that the inertia force, damping force, and stiffness force are nonlinear. However, the order of magnitude term $(N3\rho A+N42Mm)x2$ in the inertia force is relatively small. As such, the $(N3\rho A+N42Mm)x2$ term is neglected in the following sections.

## 3 Displacement-Measurement Restoring Force Surface Method

The restoring force surface method is employed to plot the internal restoring force of a nonlinear element as a function of displacement and velocity amplitudes. The nonlinear stiffness and damping forces can be visualized by cutting the restoring force surface with a vertical plane where either the velocity or displacement is near zero, respectively. Theoretically, it is better to measure displacement, velocity, and acceleration at each sampling point simultaneously. In practice, measuring the acceleration and then using an integration procedure to obtain the other two variables is common.

*f*

_{k}(

*x*).

*p*, so the tip displacement could be calculated using the modal shape function. Then, taking the three-point numerical differentiation program given in Eq. (18), the velocity and acceleration can finally be obtained as

*u*denotes the location of the sampling point, and Δ

*t*is the time interval.

It should be noted that the response trajectory of the oscillators needs to cover the phase space as much as possible to provide both enough and continuous data sets [24,25]. Therefore, the choice for the excitation source should be careful, especially for multistable structures. Using the above procedure, the restoring force surface can be defined by the triplets $(x,x\u02d9,fnon(x,x\u02d9))$. Generally, cross sections of the restoring force surface along the axes where the velocity and displacement are equal to zero yield nonlinear restoring force curves and damping force curves, respectively.

However, if the sampling data around the zero velocity are not sufficient for characterization, it will lead to a failure or cause a loss of identification accuracy. Therefore, a superior approach is presented, and the detailed process of its mathematical explanation is as follows.

*m*and

*n*are the order of displacement

*x*and velocity $x\u02d9$, respectively, and the coefficients

*α*

_{ij}can be obtained using the least-squares method. An essential procedure in the displacement-measurement restoring force surface method is cutting the restoring force surface with two vertical planes around the zero velocity instead of using just one zero-velocity plane. It is indicated that the velocity range is chosen as $|x\u02d9|\u2264\delta 1$, where

*δ*

_{1}represents the distance between the two vertical planes. Finally, by projecting all the data sets $(x,x\u02d9,fnon(x,x\u02d9))$ to the zero-velocity plane, and the nonlinear restoring force can be expressed as

The *k*_{i} can be estimated by the least square method.

*x*| ≤

*δ*

_{2}. Finally, projecting all the data sets $(x,x\u02d9,fnon(x,x\u02d9))$ to the zero-displacement plane, and the linear damping force can be expressed as

Least-squares parameter estimation can be utilized to obtain the values of the damping coefficients *c*.

*f*

_{k}(

*x*) and the theoretical one, the normalized mean-square error (NMSE) indicator [27] is defined as

*N*is the total number of samples.

The flowchart for identifying strongly nonlinear structures based on the displacement-measurement restoring force surface method is shown in Fig. 2.

## 4 Numerical Investigations

In this section, quasi-zero stiffness and bistable stiffness structures are utilized to identify the displacement-measurement restoring force surface method. Then, a comparison between the displacement-measurement and acceleration-measurement restoring force surface is discussed. Finally, the influence of the noise level on the nonlinear restoring force identification accuracy is investigated.

The selection and calculation of the fundamental physical parameters of the beam are shown in Table 1. In the following simulations, the quasi-zero stiffness and bistable stiffness beam with the function of *f _{non}* = 10

*x*+ 10

^{5}

*x*

^{3}+10

^{8}

*x*

^{5}and

*f*= −30

_{non}*x*+ 3 × 10

^{5}

*x*

^{3}− 10

^{8}

*x*

^{5}are chosen.

Symbol | Numerical values | Symbol | Numerical values | Symbol | Numerical values |
---|---|---|---|---|---|

E | 1.98 × 10^{11} N/m^{2} | M_{m} | 4.2 g | N_{4} | 8.2247 |

L | 0.15 m | c | 1 N/(m/s) | N_{5} | 901.936 |

b | 0.01 m | N_{1} | 0.034 | N_{6} | 24,727 |

h | 0.0008 m | N_{2} | 0.0545 | N_{7} | 1.356 × 10^{6} |

ρ | 7810 kg/m^{3} | N_{3} | 1.8401 |

Symbol | Numerical values | Symbol | Numerical values | Symbol | Numerical values |
---|---|---|---|---|---|

E | 1.98 × 10^{11} N/m^{2} | M_{m} | 4.2 g | N_{4} | 8.2247 |

L | 0.15 m | c | 1 N/(m/s) | N_{5} | 901.936 |

b | 0.01 m | N_{1} | 0.034 | N_{6} | 24,727 |

h | 0.0008 m | N_{2} | 0.0545 | N_{7} | 1.356 × 10^{6} |

ρ | 7810 kg/m^{3} | N_{3} | 1.8401 |

### 4.1 Comparison Between Differentiation and Integration.

An 80-second forward-sine sweep followed by a backward sine sweep signal is performed to excite the quasi-zero stiffness system. The frequency-swept range is 0–20 Hz, and the sampling frequency is 500 Hz. The classic fourth-order Runge-Kutta is adopted to solve the second-order differential equation.

Figure 3(a) shows the complete time-domain displacement response of the oscillator. Next, a 50 Hz Butterworth low-pass filter and a three-point numerical differentiation algorithm are introduced. After one and two times of conducting this process, both the velocity and acceleration responses are obtained. The resulting data sets for displacement, velocity, and acceleration are substituted into Eq. (17), and then the restoring force surface is constructed as depicted in Fig. 3(b). Meanwhile, two vertical zero-velocity planes with *δ*_{1} = 0.02 are selected to cutting the restoring force surface to obtain the intercepted nonlinear restoring force trajectories. Similarly, the intercepted damping force trajectories can be obtained by two vertical zero-displacement planes with *δ*_{2} = 0.0002. Figure 3(c) shows the comparison between the identified nonlinear restoring force and the theoretical one. There is an excellent fit with the NMSE of 0.041%. It can be viewed from Fig. 3(d) that the identified linear damping force trajectories do not match as well with the theoretical one due to the assumed small-displacement cutting planes causing upper and lower deviation of the intercepted nonlinear damping force data sets. The relative error of this case is 8.26%. An integration procedure from the acceleration time-domain sequence is also carried out for the same quasi-zero stiffness oscillator with the same excitation source. Figure 3(e) presents the acceleration, velocity and acceleration responces, and the trapezium integration and the second-order trend-term removal algorithms are adopted here to obtain the velocity and displacement. The acceleration-measurement restoring force surface shown in Fig. 3(f) is constructed to extract the nonlinear restoring force and damping force trajectories. It can be observed in Fig. 3(g) that the identified nonlinear restoring force error is only 0.087%. Nevertheless, the identified damping force has a relative error of 11.93%, as shown in Fig. 3(h).

The above results indicate that for monostable quasi-zero stiffness structures, both the displacement-measurement and acceleration-measurement methods for constructing the restoring force surface could identify the nonlinear restoring force well. Both methods need to be used with caution for damping coefficient identification because of the relatively higher error. There are nonlinear damping components in the dynamic equations derived from the magnetically coupled cantilever beams. The simulated displacement signal used in this paper is generated by this nonlinear equation. However, the identified damping force here is simplified to linear damping force. So, in Sec. 4.2, the identification of damping is not considered.

Bistable structures can exhibit both inter-well or intra-well motions dependent on the excitation and the initial conditions. Based on the contributions of Worden [24,25], the excitation source must make the response trajectory cover the phase plane as much as possible. If the bistable structure oscillates completely across the potential well or only makes intra-well motions, it is impossible to construct an excellent restoring force surface, and it will make the fitting accuracy poor. In this example, a forward-sine sweep followed by a backward sweep with a frequency range of 5–25 Hz is adopted. The swept time is 80s with a sampling frequency 500 Hz. Figure 4(a) shows the displacement, velocity, and acceleration responses of the oscillator. A 5–50 Hz Butterworth bandpass filtering program and differential procedure are conducted to obtain the velocity and acceleration sequences. In Fig. 4(b), the restoring force surface is constructed and the intercepted nonlinear restoring force and damping force are marked by circle lines. It can be observed from Fig. 4(c) that the identified nonlinear restoring force fits well with the theoretical curve, and the NMSE error is only 0.093%. However, the identified damping force shown in Fig. 4(d) still has low accuracy. As mentioned earlier, the displacement-measurement restoring force surface method can identify the nonlinear restoring force well in bistable structures, especially the nonlinear stiffness force. In Fig. 4(e), the acceleration response sequence is first presented and then make the integration procedure. Finally, the velocity and displacement sequence are given. It can be observed from Fig. 4(f) that the displacement to velocity integration procedure is reasonable, but the integrated displacement based on velocity is wrong, as shown in Fig. 4(g). Many trend-term removal methods such as mean removal, second-order trend removal, and a high-pass filter fail to obtain a reasonable displacement response. It is not easy to make the integration and detrending program to locate the inter-well and intra-well motion in the dynamic response of bistable oscillators. It can be observed in Fig. 4(h) that the constructed restoring force surface data sets are wrong.

From the aforementioned numerical analysis, it is found that the displacement-measurement restoring force surface method may be more desirable for the identification of quasi-zero and bistable structures. Therefore, it is necessary to start from the time series of the measured displacement to formulate identification procedures for bistable structures. However, the differentiation calculation of the measuring displacement will lead to noise amplification in the velocity and acceleration. In Sec. 4.2, the influence of the noise level on the identification accuracy of the nonlinear restoring force is discussed. The damping force is not discussed as no noise pollution in the identification process has led to more significant errors.

### 4.2 Influence of the Noise Level on Identification Accuracy.

In this section, the influence of the noise level on the identification accuracy of the nonlinear restoring force in quasi-zero stiffness and bistable structures is discussed. Noise with intensities of 40 dB, 30 dB, and 20 dB are added to the displacement response, respectively, and then, comparative discussions are performed.

In Fig. 5(a), the frequency-swept displacement response with 20 dB noise of the quasi-zero stiffness structure is shown, and then, the same procedure as in Sec. 4.1 is adopted for obtaining the velocity and acceleration. Figure 5(b) presents the constructed restoring force surface and intercepted nonlinear restoring force data sets. The identified nonlinear restoring forces with the noise levels of 40 dB, 30 dB, and 20 dB are shown In Fig. 5(c), and the corresponding NMSE errors are 0.23%, 1.6%, and 11.3%, respectively. The identification accuracy decreases with the noise-level increase, and identification results are close to the actual value when the noise level is below 30 dB. The 40 dB, 30 dB, and 20 dB noise levels are also added to the frequency-swept displacement response for bistable structures. Figure 5(d) presents the displacement response with 20-dB noise, and the resulting velocity and acceleration obtained by the differentiation procedure are also given. In Fig. 5(e), the constructed restoring force surface for the bistable oscillator and the nonlinear restoring force trajectory is shown. The final identification results are plotted in Fig. 5(f). Similarly, the identified nonlinear restoring force can keep a good fit when the noise level is lower than 30 dB. In general, the identification accuracy of the bistable oscillators is lower than that of the quasi-zero stiffness structures. Considering that the inter-well and intra-well chaos motion is more complicated than the frequency-swept response in monostable structures, the lower accuracy in identifying bistable structure is reasonable.

## 5 Experimental Validations

To verify the displacement-measurement restoring force surface method for nonlinear restoring force identification of three strongly nonlinear structures in realistic conditions, an experimental test rig is set up. The quasi-zero stiffness, bistable stiffness, and tristable stiffness structures are designed in the cantilever-beam system using rotatable magnetic couplings. The excitation source is generated by a vibration exciter (JZK-50, Econ Technologies Co., Ltd., Hangzhou, China) and a power amplifier (YE5874A, Econ Technologies Co., Ltd., Hangzhou, China); a vibration controller (VT-9002-1, Econ Technologies Co., Ltd., Hangzhou, China) is used to realize the feedback and real-time control of the excitation signal. A displacement sensor (HL-G112-A-C5, Keyence, Osaka, Japan) is applied to measure the relative displacement response, and the data sets are collected by an oscilloscope (MSOX3052A, Agilent, Beijing, China). Another acceleration sensor (EV4540, Econ Technologies Co., Ltd., Hangzhou, China) is installed on the fixture to monitor the vibration exciter output signal and transmit it to the data acquisition instrument (MI-70080, Econ Technologies Co., Ltd., Hangzhou, China).

The cantilever beam is made of 65Mn spring steel in the experiments and has dimensions of 146.5 × 10 × 0.8 mm. All the cantilever beams have two endmost magnets except for tristable beam 2 (low potential well tristable had four magnets). The magnet's size is 10 × 10 × 2.5 mm^{3}, and they are made from NdFeB. To extract the continuous nonlinear restoring force curve with the displacement-measurement restoring force surface method discussed in Sec. 3, the mode of forward-sine swept-frequency followed with the same intensity of backward swept-frequency is adopted for the experimental excitation. More detailed parameters of the configuration of the four typical beams and excitation levels are illustrated in Table 2.

Types | Equilibrium position (mm) | Frequency range (Hz) | Excitation acceleration (m/s^{2}) |
---|---|---|---|

Quasi-zero stiffness | 0 | 5–25 | 3.5 |

Bistable stiffness | −8.16; 0; 10.4 | 10–30 | 8 |

Tristable stiffness1 | −11.4; −6.29; 0; 5.9; 13.4 | 10–30 | 10.5 |

Tristable stiffness2 | −11; −6.7; 0; 8.2; 14.1 | 8–28 | 8 |

Types | Equilibrium position (mm) | Frequency range (Hz) | Excitation acceleration (m/s^{2}) |
---|---|---|---|

Quasi-zero stiffness | 0 | 5–25 | 3.5 |

Bistable stiffness | −8.16; 0; 10.4 | 10–30 | 8 |

Tristable stiffness1 | −11.4; −6.29; 0; 5.9; 13.4 | 10–30 | 10.5 |

Tristable stiffness2 | −11; −6.7; 0; 8.2; 14.1 | 8–28 | 8 |

The force gauge (M5-2, MARK-10 Corporation, Beijing, China) is installed on the structure driven by the ball screw, and the resolution is 0.002 N. As the dynamometer pushes the cantilever beam to deform, the movement displacement is recorded by a digital indicator. Each measurement needs to be conducted from the steady-state point position, and then, the force–displacement data are arranged from negative to positive deflection.

The linear stiffness of the cantilever beam is measured first, and the value is *K* = 80 N/m. Then, a swept-frequency excitation signal is conducted on the linear cantilever beam, and the resonant frequency is 18.5 Hz. Using the equation of *M* = 80/(18.5 × 2*π*)^{2}, the equivalent mass is *M* = 5.9 × 10^{−3} kg.

Generally, the dynamic response needed to be reconstructed based on the identified nonlinear restoring force and then compares with the measured data sets [34–36]. However, for strongly nonlinear structures like the bistable system investigated in this paper, it is almost impossible to do so. The dynamic response of the bistable system is very sensitive to the nonlinear coefficients, initial conditions, and excitation level. As such, slightly incorrect identification results will cause the oscillator to exhibit different forms of chaotic movement. In this paper, the identified nonlinear restoring force is compared with the directly measured one. Besides, the nonlinear restoring force of the actual magnetically coupled cantilever beam is not a perfect polynomial or a smooth spline, and it could not be accurately fitted. Owing to this, qualitative comparisons are considered in the experimental investigations.

The quasi-zero stiffness beam is excited by a forward-sine sweep (5–25 Hz) followed by a backward sine sweep (25–5 Hz) with a sampling frequency of 500 Hz at the constant acceleration of 3.5 m/s^{2}. Figure 6(a) shows the constructed restoring force surface, and the nonlinear restoring force curve is extracted by assuming the velocity $|x\u02d9|\u22640.02m/s$. The extracted nonlinear restoring force curve is not monotonically increasing, and both hardening and softening properties can be observed. It is difficult to construct a perfect zero stiffness zone around the steady-state point in real applications, and instead, a relatively close zero stiffness is achieved by coupling two negative stiffness regions. If the magnetic force keeps increasing and the nonlinear restoring force crosses the zero points, the bistable and tristable phenomenon will occur. Figure 6(b) shows that the identified nonlinear restoring force curve has a good agreement with the measured force–displacement trajectory.

The bistable beam is excited under a forward sine sweep (10–30 Hz) followed by a backward sine sweep (30–10 Hz) with a sampling frequency of 500 Hz at the constant acceleration of 8 m/s^{2}. Figure 6(c) shows the constructed restoring force surface of bistable beam and by assuming the velocity $|x\u02d9|\u22640.02m/s$, the nonlinear restoring force curve can be extracted. It can be observed from Fig. 6(d) that the identified nonlinear restoring force by the proposed method is reliable when compared to the measured force–displacement trajectory.

The potential wells and equilibrium points are critical factors for designing a multistable vibration structure or putting it into practical use [37,38]. For example, a tristable piezoelectric energy harvester needs to have a nonlinear restoring force design based on environmental excitation's intensity and frequency range. Moreover, the equilibrium points cannot be designed to be too large as this will cause piezoelectric layers to break up.

To verify the ability of the displacement-measurement restoring force surface method for nonlinear restoring force identification of different multistable structures, two asymmetric tristable beams with high and low potential wells are compared in Fig. 7. Both tristable oscillators’ restoring force surfaces are constructed well, and the nonlinear restoring force can be extracted, as shown in Figs. 7(a) and 7(c). Consequently, by comparing the identified results with the measured force–displacement trajectories in Figs. 7(b) and 7(d), the proposed identification method has strong applicability to multistable structures with different potential wells and equilibrium points.

## 6 Conclusions

This paper investigates the displacement-measurement restoring force surface method to identify the nonlinear restoring force in strongly nonlinear structures. To verify the effectiveness and accuracy of the identification procedure, quasi-zero stiffness and bistable structures are numerically investigated. By comparing it with the acceleration-measurement restoring force surface method, it is found that the method of measuring acceleration and performing the integration procedure cannot obtain a reasonable inter-well and intra-well displacement response in bistable systems. So, starting from the time series measurement of displacement to construct the restoring force surface for nonlinear restoring force identification is necessary. In the investigated monostable quasi-zero stiffness structure, the displacement-measurement method provides more accurate results if no noise is considered. Additionally, the displacement-measurement restoring force surface method is robust for identifying the nonlinear restoring force in strongly nonlinear structures when the noise level is lower than 30 dB. In experimental conditions, quasi-zero stiffness, bistable, and tristable cantilever beams are designed to obtain their displacement response for identifying the restoring force. The identified nonlinear restoring force curve has a good agreement with the directly measured force–displacement trajectory. Finally, the identified results of two asymmetric tristable structures show that the displacement-measurement restoring force surface method has strong applicability for multistable structures with different potential wells and equilibrium points.

## Funding Data

National Natural Science Foundation of China (Grant No. 51975453; Funder ID: 10.13039/501100001809).

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request. The authors attest that all data for this study are included in the paper.