Design and characteristics of a novel QZS vibration isolation system with origami-inspired corrector

A new class of quasi-zero-stiffness (QZS) vibration isolation systems, inspired by the origami metamaterial, is proposed to achieve high-performance vibration suppression in this paper. According to the mechanical characteristics of Tachi–Miura origami (TMO) with single degree-of-freedom, the nonlinear geometric relationship is developed with the folding angle as the master variable. By utilizing equivalent transformation and virtual work principle, the static model is established, the influence of structural parameters on stiffness is investigated, and the negative stiffness mechanism of origami mechanism is revealed. By adding a linear spring with positive stiffness to the origami in parallel, the Tachi–Miura origami vibration isolator (TMOriVI) is obtained. Subsequently, the governing equation is presented by means of the harmonic balance method. Two types of instability situations, jump phenomenon and unbounded response, are studied, and their analytic criteria and relationship are derived. Finally, through the parametric influence analysis and comparison with existing QZS isolator, the effectiveness and superiority of the proposed isolator are verified. The proposed vibration isolation system with great design flexibility exhibits a significant potential in the field of low-frequency vibration isolation.


Introduction
Vibration phenomena exist widely in daily life and industrial production. Although there are some examples of taking advantage of vibration [1][2][3][4], most of the vibrations are harmful. In mechanical equipment, vibrations tend to trigger serious mechanical faults [5,6]. The altering stress generated by the vibration can induce material fatigue failure for the mechanical parts and supporting structures. For precision equipment, vibration is likely to cause poor machining quantity and unsatisfactory measurement accuracy [7]. Furthermore, vibration even brings about hidden danger to human health and safety to a certain extent [8]. Therefore, achieving efficient vibration isolation has always been the pursuit in industry.
Roughly speaking, the vibration isolation techniques can be grouped into three categories: active, semi-active, and passive vibration isolation systems [9,10]. The active vibration system is comprised of the servo mechanism with sensors, the signal processor, and the actuator, which relies upon the external energy to control the vibration. The semi-active vibration isolation system is realized by the components which can adjust the stiffness coefficient or damping coefficient actively.
Because of the concise structure, easy implementation, and high-cost performance, the passive vibration isolator has been favored by the engineering circles for many years. Traditional linear passive vibration isolation systems are developed by introducing a linear spring element and a linear damping element between the isolated object and the vibration source, which enable the external energy to be stored and then dissipated. Unfortunately, affected by the inherent property of the system, this technology can only take effect when the frequency of the vibration source is √ 2 times of the isolator's natural frequency [11], leading to the failure to achieve low-frequency vibration isolation. Obviously, deducing the natural frequency will render the static deformation larger, which also needs to be avoided in practical applications. To manage this dilemma, the high-static-lowdynamic-stiffness (HSLDS) vibration isolation concept was introduced to the field in [12], namely, this kind of isolator can realize large static load and ensure small dynamic deformation, simultaneously. Particularly, when the dynamic stiffness at static equilibrium is equal to zero and other stiffness is always greater than zero, this characteristic is called quasi-zero-stiffness (QZS). Based on this idea, a large number of highperformance HSLDS vibration isolators [13][14][15][16][17] with various principles and configurations have been proposed and investigated. For this kind of nonlinear isolators, the key factor is to design the negative stiffness mechanism [18][19][20][21][22], then, which can be used to optimize or combine to form HSLDS isolator. Among these technologies, the approaches to achieve negative stiffness can be summarized as follows.
The most common form is the energy storage recovery mechanism. Euler strut is a basic device providing negative stiffness. When the axial load on both ends of the member bar exceeds the critical load, the radial stiffness of the bar is negative. Platus [23] designed a vibration isolator with vertical and horizontal stiffness close to zero based on this principle. More typically, Carrella presented the three-spring vibration isolator with symmetrical oblique springs [12,24]. The preload is applied until the oblique springs is perpendicular to the main spring, and this position is taken as the static equilibrium point. The investigation showed that this three-spring configuration can exhibit QZS characteristic and good isolation performance. Furthermore, Huang et al. combined the compressed Euler buckled beam and three-spring mechanism to propose a new type of vibration isolator, and analyzed the isolation performance under different excitations in detail [25,26]. Note that the QZS vibration isolators based on magnetic force [27,28] or pneumatic force [29,30] can also be regarded as a special kind of energy storage recovery mechanism. When they deviate away from the equilibrium position after disturbed, a force in the same direction as the displacement is generated to achieve negative stiffness.
Another effective means to reach negative stiffness is utilizing the nonlinearity generated by the structural shape or the mechanism motion. The fairly typical case is the cam mechanism. Li et al. designed and fabricated two different HSLDS vibration isolators with a cam-roller-spring mechanism [31,32] to test the isolation performance in the low-frequency range. Inspired by the bird limb skeleton, Jing proposed a series of HSLDS vibration isolation systems based on nonlinear scissor-like structured platforms and combined the parallel mechanism to expand the vibration signal to multidimensional vibration isolation platforms [33][34][35][36]. Moreover, the designed isolation systems were tested in different application scenarios, such as post-capture mechanism [37], hand-held jackhammer [8], vehicle suspension [38], and tracked mobile robot [39], and delivered outstanding isolation performance.
Besides, an emerging methodology is to adopt metamaterial, a class of materials or structures with counterintuitive mechanical properties. Trishan [40] designed the metamaterial comprised of regular array of interlocked rigid hexagonal subunits, which could display simultaneous negative Poisson's ratio and negative stiffness. Chen [41] proposed the lattice metamaterials with negative Poisson's ratio by using sinusoidally curved beams. This characteristic endows the metamaterial with remarkable broadband vibration mitigation capability. Tan [42] designed a cylindrical negative stiffness metamaterial composed of spatial curved beams, and the impact test results exhibited good cushion performance.
On the other hand, origami, an ancient folk handicraft art, has received extensive research attentions with attractive characteristics. One point should be emphasized, that is, origami mechanism is actually a special kind of metamaterial with negative Poisson's ratio and negative stiffness in a certain range [43].
Benefited from the abundant folding deformation and the corresponding unique mechanical properties, origami has gradually developed into a wide potential engineering technique in the areas of space antennas [44], mobile robots [45], actuators, developed mechanisms, energy-absorbing devices [46], etc. However, limited by the complicated mathematical model and extremely strong nonlinearity corresponding to the special configuration, origami is still in the initial stage. The proposed designs just consider the kinematic characteristics of origami, not the dynamics. And there are few examples of using it as the isolator directly. The aforementioned obstacles motivate us to propose a new way to realize effective vibration suppression by employing origami as the vibration isolator.
Inspired by the origami metamaterial, a novel passive vibration isolator with an origami corrector is proposed in this paper. The negative stiffness in this system is generated by the oblique creases. By tuning the geometric parameters of the origami mechanism to gain an appropriate negative stiffness and paralleling a positive stiffness linear spring as the compensation to keep the system stable, then, the HSLDS or even QZS can be achieved. In addition, the stability of the dynamic response considered here concludes two situations, specifically, unbounded response and jumping phenomenon. The analytic criteria are derived based on the property of the amplitude-frequency equation, which can accurately determine the system parameters within the allowable ranges. Furthermore, in this study, the vibration isolation performance with various parameters is evaluated, and the comparative studies are conducted to verify the advantages of the proposed isolator.
The remainder of the paper is organized as follows. Section 2 presents the origami-inspired vibration isolation system and its geometric relationships. The static modeling and negative stiffness mechanism of the proposed origami isolator are developed in detail in Sect. 3. Based on the equation of motion, Sect. 4 solves the frequency response of the isolator and investigates the stability criteria. The vibration isolation performance with various system parameters is estimated and the comparison and discussion are proposed in Sect. 5. Finally, Sect. 6 gives the conclusion.

Description of Tachi-Miura origami
As the key issue of origami, crease design directly determines the geometric shape, mechanical properties, and functional characteristics of the final product. Among origami science, the typical crease design samples include Miura, waterbomb, Yoshimura, and diagonal crease designs [45]. The origami mechanism and derivative mechanism designed according to these creases and their combination are numerous, and Tachi-Miura origami (TMO) [47] is a classic form of three-dimensional origami. Compared with other origami mechanisms, TMO has the following important characteristics. First, TMO is the combination of Miura origami and three-dimensional mechanism, and its basic crease is Miura-ori, the most widely studied and frequently used folding form. The TMO retains the geometric properties of Miura-ori and therefore shows appropriate representativeness. Second, TMO belongs to the tubular folding configuration and has better structural stability to adapt to greater bearing capacity, which is of great benefit to practical engineering. Third, the TMO is a kind of action origami and has one degree-of-freedom in the vertical direction. There is no need to realize the folding and unfolding functions through the deformation of the facets. In view of these advantages, we adopt the TMO as the origami-inspired vibration isolation corrector in this paper, and the model is shown in Fig. 1. The Miura origami is designed as the prototype to be improved, and then the 3D tubular mechanism shown in Fig. 1a is developed. Figure 1b is the 2D structure, and when several planar structures are folded by Miura-ori according to certain specifications, the three-dimensional tubular mechanism can be formed. Figure 1c is a relatively simple two-layer origami structure. Obviously, the TMO corrector is composed of several basic origami modules stacked in series in the vertical direction. Due to the symmetry of the two origami patterns, the difference between the mountain and the valley creases is only reflected in the structure of TMO, but does not affect the mechanism analysis. For the sake of simplicity, we take the upper subfigure for example to study. In Fig. 2 While folding along the crease shown in Fig. 2, the space shape illustrated in Fig. 3 can be formed. The two subfigures are the folding states of the same origami cell in different perspectives. The local coordinate frame Fig. 3. z-axis coincides with the height direction, y-axis is collinear with the middle mountain crease M 1 M 2 , and x-axis is determined by the right-handed rule. θ is defined as the folding angle, so that U 1 M 1 D 1 = 2θ , and the folding range is θ ∈ [0, π/2]. The plane angle of the dihedral corner made up of two oblique folding facets is represented as 2γ , and the folding interval is γ ∈ [0, π/2]. β is the angle between the crease M 2 M 3 and y-axis positive direction.
It is supposed that the plane Σ is developed by the creases M 2 D 2 and M 2 U 2 . Draw the vertical lines of crease M 2 M 3 through the points D 2 and U 2 , and point E is the intersection. Subsequently, make the perpendicular bisector of segment U 2 D 2 through the point E, and the foot point is marked as T , so that point T is the middle point of segment U 2 D 2 . According to the solid geometry knowledge, we have Similarly, make the normal line of plane Σ through point D 1 , and the intersection is labeled as point J . The segment I J is perpendicular to the crease M 2 D 2 , then we can obtain We describe Eqs. (2) and (3) in the angle form as follows Squaring on both sides of Eq. (5) leads to Utilizing the square relation of trigonometric function, And Eq. (5) can be rewritten as follows sin 2 γ = 1 − tan 2 α cos 2 θ sin 2 α 1 + tan 2 α cos 2 θ By dividing the two equations (7) and (8), the angle relationship can be characterized as Furthermore, by eliminating the variable β, we can find tan γ = cos α tan θ Here, as for TMO, the four basic geometric angle formulas are established as follows Besides, considering N -layer origami structure, the overall height H at any position is determined as

Interference-free condition
During the motion process of origami-type mechanism, the expansion and contraction in the height direction may cause the physical collision between the internal structures. In order to avoid the interference, we select the basic origami pattern to analyze, and the fully folded effect is shown in Fig. 4. When fully folded along the crease to θ = 0, the points U 1 , M 1 , and M 4 are transformed to the points U 1 , M 1 , and M 4 , respectively. Drawing a line perpendicular to the segment U 2 U 3 through the point U 1 , we can obtain the length of U 1 U 4 (13) Then, the interference-free condition of TMO vibration corrector is obtained as After simplification, we get the following concise expression as the interference-free condition

Static analysis
In this section, the mechanical analysis model of TMO is established by employing the virtual work principle. Subsequently, the influence of structural parameters on the restoring force and stiffness is studied in detail and the HSLDS characteristic is obtained. Furthermore, the negative stiffness mechanism of TMO is discovered by analyzing different types of creases, and some design strategies for obtaining negative stiffness are briefly outlined.

Mechanical modeling
In this subsection, combined with the kinematic characteristics of origami, some necessary assumptions and equivalent principles are given firstly. Assume that the facets will not deform in the process of origami motion, and all of the motion relations will be generated at the creases. Analogous to the ordinary kinematic pairs, the crease can be regarded as a rotating joint, and the facet surrounded by the creases is similar to the connecting rod. Considering the resistance of origami mechanism during folding or unfolding, a spring (similar to torsion spring) is introduced to simulate the antagonism. Furthermore, for convenience, the damping effect on the origami is simplified by using a linear damper.
Observing the origami mechanism shown in Fig. 1, we can divide the folds into two categories, horizontal crease and oblique crease. The former is always parallel to the x y plane in motion, while the latter is inclined to the x y plane. Therefore, the equivalent mechanical model, as shown in Fig. 5, can be obtained. In this figure, F is the driving force, and when the origami mechanism is in static equilibrium, F is also the restoring force. θ and γ are the folding angles of horizontal and oblique creases, respectively.
Suppose that the torque at the crease varies linearly with the folding angle, the crease torque magnitude generated by the horizontal crease of unit length is described as where k hc represents the torsional stiffness per unit length of horizontal crease. It is worth noticing that the origami mechanism discussed in this paper is based on rigid origami and the creases are regarded as the torsion springs with certain stiffness to meet the load-bearing requirements. As for the selection of origami materials, the key principle is that the stiffness of the facet is far greater than that of the crease, that is to say, rigid materials can be selected to simulate the plate surface, and soft materials can mimic the crease to realize the final mechanism. On the other hand, from the perspective of the vibration isolator, the vibration amplitude at the equilibrium position is relatively small, then the crease stiffness can be approximately modeled as linear regardless of the origami materials. Therefore, the effect of material on the performance can be transformed into the influence of linear stiffness at the crease. Similarly, the crease torque magnitude created by the oblique crease of unit length is expressed as where k oc denotes the torsional stiffness per unit length of oblique crease.
When the system is in static equilibrium state, by applying a virtual displacement, the system can produce a minimum displacement δu, and the angles at the horizontal crease and the oblique crease change δθ and δγ , respectively.
The principle of virtual work states that the total virtual work of all forces acting on the static equilibrium is zero for a group of infinitesimal virtual displacements from equilibrium. Then, we can obtain the equation as follows where L hc and L oc are the total lengths of horizontal creases and oblique creases of TMO, respectively. And they can be calculated as follows where M is the number of origami patterns in singlelayer structure of origami vibration isolator, and M = 4 hereafter.
We define the deformation of TMO at height direction as u, and combing Eq. (12), we have where H in is the initial height. By calculating the variations of Eqs. (20) and (10), we can obtain Combing the above equations, the restoring force of TMO can be written as Here, we should introduce some dimensionless parameters as followŝ whereû is the non-dimension displacement, denoting the ratio of folding,F is non-dimension force, η represents the ratio of crease-length, and λ is the torsional stiffness ratio of two types of creases, respectively. Meanwhile, we unite Eq. (20) and Eq. (24), and we define ν sin θ , κ cos θ , then the non-dimension displacement can be rewritten aŝ where ν in stands for the sine of the initial folding angle. Substituting Eq. (25) into Eq. (23), we can obtain the non-dimension restoring force as followŝ where Decomposing Eq. (26) into horizontal force and oblique force, we find Taking dimensional displacement derivative ofû, we can get the non-dimension stiffness caused by horizontal crease aŝ Moreover, the derivatives ofF oc with respect toû iŝ Thus, the total dimensionless stiffness of TMO can be composed aŝ

Effects of geometric parameters on static characteristics
In this subsection, we choose N = 7, α = 60 • , η = 1, λ = 1, and θ in = 80 • as the original parameters to study their influence on the static characteristics. The TMO's restoring force and static stiffness for the number of basic module stacking layers of N = 6, 7, 8, and 9 are plotted in Fig. 6. It can be seen from the restoring force results in Fig. 6a that with the layer increasing, the restoring force decreases. As shown in Fig. 6b, the absolute value of stiffness is inversely proportional to the number of layers. Combing Eqs. (27) and (30) leads tô The result is also consistent with the fact that the TMO can be equivalent to N -layer origami basic modules in series in the vertical direction.
Besides, when other parameters are fixed, the stiffness varies from positive to zero to negative, and in the shadow zone in Fig. 6b, the TMO exhibits negative stiffness characteristic, furthermore, which is independent of the layer number. Figure 7 displays the static analysis results for crease angles of α = 65 • , 60 • , 55 • and 40 • . It can be observed that with the increase of crease angle, the TMO is gradually transformed from a pure positive stiffness mechanism to a multi-stable mechanism with negative stiffness. In addition, the crease angle is larger, the initial positive stiffness value is bigger, and the negative stiffness characteristic is more apparent. Different from the layer number, the angle value directly affects the interval of negative stiffness.
The characteristic curves of TMO with various crease stiffness ratios and total length ratios are shown in Figs. 8 and 9. The effects of both on the TMO static characteristic are the same, which can be verified according Eq. (26). When the values of the two ratios are adjusted, the negative stiffness range will change accordingly with other parameters fixed. And the larger the ratio value, the larger the negative stiffness range. As the stiffness ratio or length ratio increases, the absolute value of stiffness increases, which is the exact opposite of the layer number.
The influence of different initial folding angles θ in on the dimensionless restoring force characteristics and dimensionless stiffness characteristics of the mechanism is shown in Fig. 10. As the initial angle increases, the static nonlinearity of TMO become grad- ually stronger. Once the angle exceeds a certain threshold, the TMO presents negative stiffness characteristic, and the stiffness nonlinearity is directly related to the crease angle value. Compared with other parameters, it can be found that TMO is more sensitive to the initial folding angle. This phenomenon can be explained by the result that the folding angle is the direct factor affecting the stiffness, but other parameters should change the folding angle to indirectly affect the stiffness.

Negative stiffness mechanism
In this part, the negative stiffness mechanism of TMO is analyzed. According to Eqs. (27) and (30), the stiffness can be decomposed into two parts, i.e., horizontal stiffness and oblique stiffness. Transforming horizontal stiffness Eq.
whereK * hc = cos θ + (θ − θ in ) sin θ . Since the coefficient 4 sin θ in /N cos 3 θ is always positive, the sign of horizontal stiffness depends on K hc . By differentiating it, we havė It is easy to find that if θ = θ in , thenK * hc attains a minimum value cos θ in , which is always positive on the domain. Thus, the horizontal stiffness of TMO is consistently positive.
By transforming oblique stiffness expression Eq. (29), we obtain K oc = 4ηλ sin θ in cos α N cos 3 θ cos 2 α + sin 2 α cos 2 θ 2K * oc (34) whereK * oc = (γ − γ in ) cos 2 α sin θ + 3 sin 2 α cos 2 θ sin θ ) + cos α cos θ , and the stiffness sign determined byK * oc . When the oblique stiffness is negative, we get the further simplification 1 + (γ − γ in ) 1 + 3 cos 2 θ tan 2 α tan γ < 0 This is the necessary condition of exhibiting negative stiffness for TMO mechanism. Two examples are given to illustrate the effectiveness of the obtained negative stiffness mechanism, and the results are plotted in Fig. 11. We only change α, then the horizontal forces are the same, and the oblique forces are different.
To sum up, for TMO, the horizontal crease plays the role of load bearing and does not possess negative stiffness characteristic. Through reasonable design of structural parameters, the oblique crease can achieve negative stiffness when the condition of Eq. (35) is satisfied.
In view of this, we can obtain two methods to design high-performance vibration isolator. Firstly, by optimizing the structural parameters, the negative stiffness caused by the oblique crease can be eliminated, so that the whole TMO has quasi-zero-stiffness characteristics, so as to achieve the purpose of high-performance vibration isolation. Secondly, we can parallel a spring with positive stiffness on the original TMO, i.e., we make full use of the negative stiffness of the origami mechanism to achieve vibration isolation in a larger range. The latter regards the origami mechanism as a stiffness corrector, which is more simple and ingenious in engineering applications. Moreover, it is worth noticing that the results of negative stiffness are consistent with the conclusions in [47,48] 4 Dynamic modeling and analysis Figure 12 displays the dynamic model of Tachi-Miura origami vibration isolator (TMOriVI) under base excitation. Suppose that the isolation system is under the base excitation of w, and the absolute displacement of the isolated object is z, then the related displacement can be represented as u = z − w. The origami isolator can be regarded as the combination of equivalent spring and equivalent damping. In the paper, for the sake of simplicity, the damping is treated as a linear element.

Equation of motion
Ignoring the mass of vibration, and selecting the absolute displacement of isolated object as the generalized coordinate, we can describe the overall system's kinetic energy as where m is the sprung mass. Now consider the static displacement of TMOriVI. The overall system's potential energy should include three terms: the elastic potential energy generated by static displacement, the elastic potential energy generated by the parallel spring, and the gravitational potential energy of sprung mass. The result can be written as where k ad is the stiffness coefficient of the parallel spring, and g is the gravitational acceleration. The Rayleigh's dissipation function can be obtained as where c is the coefficient of viscous damping. Utilizing the Lagrange equation, we have d dt where L = T − V is Lagrangian function. Substituting Eqs. (36), (37) and (38) into Eq. (39) to obtain Since the gravity of isolated object is equal to the elastic potential energy generated by the static displacement, the following equation holds.
And the equation can be rewritten as Introducing the variableū = u − u st as the equilibrium displacement, we can rewrite Eq. (43) as mü + cu + k adū + 4L hc k hc (Θ + Ξ ) Recalling the geometric relationship of TMO, the height of TMOriVI can be described as We describe Eq. (44) as the Taylor expansion of the equilibrium displacementū at the equilibrium configuration as follows where the coefficients are given in Appendix A. Substituting the coefficients into Eq. (44), we find Thus, the following equivalent form can be obtained Since the stiffness coefficient of the parallel spring is equal to the absolute value of the maximum negative stiffness, the quadratic coefficient must be zero, and the proof is detailed in Appendix B.
Therefore, the governing equation of TMOriVI under base excitation can be established as According to the above analysis, the approximate restoring force curve and the theoretical restoring force curve of TMOriVI are drawn in Fig. 13. It can be seen from the figure that at the static equilibrium, the approximate curve via Taylor expansion can be well consistent with the theoretical one. With the increase of displacement amplitude, the error between the two gradually increases, and the coincidence effect is acceptable within the displacement range of ±0.05, which also satisfies the actual engineering requirements.

The steady-state response
In this part, we use harmonic balance method (HBM), an efficient way of dealing with nonlinear problems [15,49], to obtain the steady-state solution .
Suppose that the base excitation is harmonic, that is, Introduce the following non-dimension variableŝ where (·) denotes the non-dimension time τ derivative. By employing HBM, we can develop the fundamental solution into the Fourier series with the frequency of Ω, retain the first-order harmonic terms, and obtain the solution of Eq. (55) as followŝ Introducing an intermediate variable φ = Ωτ + ϕ to rewrite Eq. (55) as Omitting the higher-order harmonic terms and letting the coefficient of the first-order term be equal to zero, we have Eliminating the variable ϕ to derive the amplitudefrequency relationship of TMOriVI under base excitation as Eliminating the variable W from Eq. (58), we can obtain the phase-frequency relationship as follows tan ϕ = − 8ξΩ 3ρ 3 A 2 + 4ρ 1 − 4Ω 2 (60)

The stability of the steady-state solution
It can be found from Eq. (55) that the governing equation of TMOriVI is consistent with that of Duffing oscillator. However, Duffing oscillator received external excitations perhaps generates some unstable anomalies, including jump phenomenon and unbounded response. The former will jump in the process of frequency sweeping test, worsening the vibration isolation effect and even damaging the isolated equipment. The latter will generate infinite displacement in the process of vibration, leading to poor isolation performance. Therefore, both of them are fully considered and avoided in the design process. For Duffing oscillator, the jump point locates at the point with infinity slope, that is, If Eq. (62) is unsolvable in the domain, then the dynamic system is stable with no-jump. If Eq. (62) has only one solution, i.e., the solutions in the domain are equal, then the system is in a critical jump state. And if Eq. (62) has two different feasible solutions in the domain, then there is a jump phenomenon in the dynamic system. Then, we can infer further that, for the critical jump state, the second-order derivative is equal to zero at the critical point.
According to the inference, taking amplitude derivative of frequency, we can get We continue to take the second derivative to obtain For the jumping points, letting Eq. (63) be equal to zero, we have By uniting Eq. (66) and Eq. (65), the relationship between amplitude and frequency at the critical state can be solved as follows The result hints that the determine of critical point is only related to the stiffness, but independent of damping and external excitation.
By substituting Eq. (67) Replacing the non-dimensional variables in Eq. (71) with original expressions in Eq. (54), we can get the necessary and sufficient condition of critical jump state for TMOriVI as follows Thus, the jump criteria can be summarized as follows: (1) If L P > 0, then there are two jump points for TMOriVI; (2) If L P = 0, then the dynamic system is in a critical state and is essentially stable; (3) If L P < 0, then there is no jump phenomenon for TMOriVI.
Once the proposed vibration isolator is designed, the dynamic parameters including k 1 , k 3 , m, and c are all determined, then the maximum excitation amplitude in the stability range can be calculated as Therefore, when the base excitation amplitude is no larger than W M e , the dynamic system is stable. In particular, if TMOriVI has QZS feature, then k 1 = 0, we can get the following concise expression When the amplitude of base excitation W e is greater than W Z e , there is the jump phenomenon in TMOriVI dynamic system, and vice versa.
The above criteria give the boundary jump-avoidance boundary, which can be used in the design process to ensure a stable and efficient vibration isolation system. In addition, the vibration isolation system subjected to base excitation may has unbounded response, a unstable response case that the peak value of amplitude response is infinite. When the unbounded response occurs, the peak amplitude only matches one frequency, then we have If the unbounded response occurs, then the two real solutions of the above equations both are negative. By using Vieta theorem, the following relations can be obtained Solving the above equations, we get the critical condition of unbounded response as Once the excitation amplitude is greater than W ub , there is the unbounded response in the dynamic system. It can be found that the decisive factors of unbounded response only include ρ 3 , ξ and W , yet are independent of ρ 1 .
And substituting Eq. (78) into Eq. (71) to find Similarly, the maximum amplitude of base excitation for the dynamic system without unbounded response can be calculated as The above results turn out that unbounded response is a special and extreme class of jump phenomena. Figure 14 shows the unstable response. We select the amplitude of base excitation as the variable, the stable, critical, jumping, and unbounded cases are plotted in this figure, and the shadow area obtained by Eq. (72) means the jumping range. With the increasing of W e , the instability of amplitude-frequency response becomes aggravated. And obviously, the unbounded response is included in the jump area, which validates the relation described in Eq. (79). the vibration isolation performance of the proposed TMOriVI. Subsequently, the effects of the vibration system parameters on the isolation capability are analyzed. Finally, the fair comparison with existing QZS isolator is presented to show the advantages of the proposed origami-type isolation system.

Evaluation indexes of vibration isolation performance
The function of the imposed vibration isolator is to suppress or dissipate the external vibration energy to protect the equipment and to improve system performance. In essence, the goal is to make the vibration from the base transmitted to the sprung mass as small as possible. In order to characterize this effect, we adopt displacement transmissibility, a widely used and effective test index, to measure the vibration isolation effect. The absolute displacement transmissibility T a d can be computed in the decibel form as follows where Am[ẑ] and Am[ŵ] denote the amplitudes of the displacement transmitted to the isolated equipment and the original harmonic displacement, respectively, and the initial phase angle ϕ can be solved according to Eq. (58). The effective initial isolation frequency is another index used extensively to assess the vibration isolation effect, and a lower effective initial isolation frequency means a wider isolation band.
Owing to jump phenomenon elaborated above, there perhaps exist piecewise amplitude, affecting the efficient initial isolation frequency directly. Considering the stability of the system, we can defined the effective initial isolation frequency as follows where Ω jm is the greater value of the two jumping frequencies, and Ω da denotes the frequency value corresponding to T a d = 0. Letting T a d = 0 leads to A + 2W cos ϕ = 0 (83) Combing Eq. (58) to get the transform relationship as follows Substituting Eq. (84) into Eq. (59), we can obtain

Parametric influence and discussion
In this subsection, the effect of parameters on the isolation performance and some discussion are presented. The default parameter settings are W = 0.04, ξ = 0.04, N = 7, η = 1, λ = 1, θ in = 80 • and α = 60 • , and the red solid curves represent the default parameters in the following figures. The effects of excitation amplitude W and damping coefficient ξ on vibration isolation performance are shown in Figs. 15 and 16. It can be seen that increasing the amplitude or decreasing the damping will increase the initial frequency of vibration isolation, which perhaps causes the curves of displacement transmissibility to bend, and further deteriorates the isolation performance. Noticing the fact that in the cases of small amplitude and large damping, the system can obtain pretty ideal isolation performance with almost full band range, as shown by the orange dashed line and blue dotted line in the figure. However, from the perspective of absolute displacement transmissibility, there is a  marked difference between the two factors. In the effective vibration isolation bandwidth, reducing the excitation amplitude can directly reduce the transmissibility without curves crossing, and effectively improves the vibration isolation performance for any same excitation frequency. It can be explained that, for the nonlinear dynamic system, the excitation amplitude directly or linearly affects the vibration isolation performance. On the contrary, for distinct damping coefficients, there is an intersection region greater than the initial isolation frequency in the transmissibility figure, rendering the curve with smaller damping lower. It means that increasing damping is not completely conducive to the improvement of vibration isolation performance, especially in the high-frequency range, which is a major difference with W . Since the crease length ratio and stiffness ratio play the completely same role in the governing equation, we simplify them as one variable termed μ = ηλ. Figures 17 and 18 exhibit the influence of the layer number N and the equivalent ratio μ, respectively. We can clearly find that the changing trends of layer number and equivalent ratio are substantially similar. As the layer number increases or the equivalent ratio decreases, the superior vibration isolation performance can be achieved. Specifically, the effective initial isolation frequency on the abscissa axis moves towards left to widen the effective isolation frequency band, and the transmissibility curve moves toward down with no bending and less absolute displacement. Similar to the excitation amplitude W , there is no intersection between various curves with different values, showing that the isolation effect corresponding to the same frequency is more beneficial. The performance results can be interpreted by Eq. (51). There is a negative cubic relationship between N and k 3 , resulting in the flattening of the restoring force cure and achieving more efficient HSLDS. The relation between μ and k 3 is approximately linear. Like layer number N , with the equivalent ratio decreasing, the restoring force curve varies more gently, thus, the better isolation performance can be reached.
The effects of two angles, initial folding angle θ in crease angle α, are presented in Figs. 19 and 20, respectively. The apparent feature is that the vibration isolation performance is negatively related to the values of the two angles, no matter in effective initial frequency or in absolute displacement transmissibility. This phe- With the decrease of θ in or α, the initial height of TMOriVI reduces synchronously, leading to smaller static equilibrium folding angle. Combing with Eq. (51), the cubic stiffness term also becomes smaller, which is the main cause of the improved isolation performance. Besides, there is no intersection between transmissibility curves with various values. It is worth noticing, however, that when the two angles are below a critical value where the minimum stiffness of TMOriVI is equal to zero, not a negative value, then there will be a sudden change in the vibration isolation curves of the dynamic system. For example, when θ in = 72 • or α = 40 • , the isolator without positive stiffness spring in parallel, i.e., k 1 > 0, k ad = 0, achieves a poor isolation performance due to the loss of the QZS characteristic. Moreover, this result is also consistent with the previous stiffness analysis as shown in Figs. 7 and 10.
Next, we consider the effect of the additional parallel spring on the vibration isolation performance. Since the case of pure positive stiffness of TMOriVI has been discussed above, herein, we are merely concerned about the isolation performance of the isolator with parallel spring. When the external positive spring exactly balances out the negative stiffness of TMOriVI labeled k * ad , namely, k 1 = 0, the vibration isolation effect reaches the optimal value, and the vibration isolation system achieves QZS characteristic. However, once the spring stiffness exceeds this critical one, such as 1.1k * ad and 1.3k * ad shown in Fig. 21, the effective initial isolation frequency and absolute transmissibility both increase, and the isolation performance becomes weakened. The principal reason for this is the increased first-order stiffness term. Because the first-order stiffness term increases, the restoring force curve becomes steep, leading to the degradation of the HSLDS characteristic. Moreover, the larger the stiffness of the parallel spring is, the more serious the performance degradation is. Meanwhile, it is worthy to emphasize that when k ad = 2.77k * ad , the effective initial isolation frequency ratio is √ 2, and the isolation performance of the origami isolator is still better than that of the linear system within the effective vibration isolation frequency band, which indeed demonstrates the effectiveness of the origami-type isolator.

Comparison with the existing QZS isolator
The most famous QZS isolator is without doubt the three-spring vibration isolator (TSVI) that has been investigated by many researchers so far [9,12,24]. In view of this, we select the TSVI as a reference. To compare the vibration isolation effect between TMOriVI and TSVI impartially, we consider the isolators with the same loading mass and the same natural frequency. For TMOriVI, the dynamic parameters in the calculation are listed in Table 1.  Figure 22 displays the model of TSVI. In this isolation system, by selecting the parameters, the equivalent stiffness in the vertical direction can be adjusted to a very small value. The restoring force in the vertical direction can be given as where L in is the initial length of the horizontal spring. The force-displacement relationship at the static equilibrium position can be simplified as follows For the system, the governing equation for harmonic base excitation is obtained as The non-dimensional form is given bŷ And the corresponding transformation factors are listed aŝ We can find the condition of QZS isolator Because the damping coefficients of the isolators can be tuned by the external dampers, we set the values of  Table 2 are taken in the following numerical calculation. The displacement transmissibilities of the two QZS isolators under base excitation are compared in Fig. 23. In general, both of the two isolators can achieve significant vibration isolation effects, and their isolation performance is almost the same at high frequency, which is due to their QZS characteristics. However, compared with TSVI, the proposed vibration isolator has some advantages. First, the proposed TMOriVI holds better stability. It can be seen from the displacement transmissibility curve that the TMOriVI is stable and has weak nonlinearity, without the multi-steady phenomenon. By comparison, the jump response occurs for the TSVI, which may bring hidden risks for practical engineering applications. Second, the effective initial isolation frequency of the TMOriVI is much lower, and a wider isolation band can be achieved. It is observed from the figure that the TMOriVI has a smaller abscissa intersection value, let alone considering the influence of jumping, which means the proposed origami-type isolator can suppress large amplitude vibration effectively in a wide frequency range. Third, the TMOriVI is superior to the TSVI in terms of vibration isolation efficiency in the low-frequency band. The curve of TMOriVI shown in Fig. 23 is always under another one, indicating the relatively more efficient vibration isolation performance.
Furthermore, we discuss the range of cubic stiffness for the two QZS isolators. For the TSVI, the condition of QZS characteristic is q 1 = 0, then the following Then, the cubic stiffness term can be rewritten as Note that the expression is strictly monotonic decreasing in its domain, and whenL = 1, q 3 has the minimum value 1. For the original expression, q 3 = 1 means the value of the cubic stiffness is equal to k v . According to the knowledge of nonlinear vibration, the smaller the cubic stiffness is, the better the low-frequency vibration isolation can be gained. However, lower vertical stiffness means the isolation system will sacrifice the loading capacity. In a word, it is an inherent defect of the system that cannot manage the tradeoff of high-load ability and ultra-low-frequency vibration isolation performance.
Moreover, for the proposed TMOriVI, the equivalent cubic stiffness can be tuned arbitrarily. For example, the other structural parameters are fixed, and we plot the nonlinear stiffness and loading mass under various single layer heights d and layer number N . The results are exhibited in Figs. 24 and 25. It can be found that the equivalent nonlinear stiffness and the loading capacity of the proposed isolation system gradually decrease with the increase of the height of the single layer. As the number of layers becomes larger, the cubic stiffness declines, while the bearing capacity can be effectively improved. By changing the values of structural parameters, the cubic stiffness of TMOriVI  can be lower than the minimum value of the TSVI, and better loading ability can be ensured. In other words, the proposed isolator can achieve superior low-frequency vibration isolation performance and greater bearing capacity at the same time.

Conclusion
Inspired by the negative stiffness characteristic of origami metamaterials, a new kind of QZS vibration isolation systems with positive stiffness spring compensation was presented in this paper. On the basis of the geometric relationship among folding angles, the folding ratio u was selected as the principal coordinate to establish the static model, and the negative stiffness mechanism was given. The equation of motion solved by HBM was used to investigate the instability modes, unbounded response and jump phenomenon, and their stability critical conditions were developed. With effective initial vibration isolation frequency and absolute displacement transmissibility as evaluation indexes, the effects of various system parameters on isolation performance were evaluated and compared with the linear passive isolator. The simulation results and comparison with the counterpart have been provided to show the effectiveness and advantages of the proposed origamitype vibration isolator. Besides, the proposed TMOriVI with many structural parameters exhibits strong designability and broad application prospect.
In the future research work, we will fabricate the vibration isolation devices to test the QZS feature and the vibration isolation performance.

Appendix B
On the one hand, according to the previous analysis of the restoring force and stiffness, the TMO belongs to multi-stable mechanism. After paralleling the spring k ad , the overall stiffness of TMOriVI is nonnegative, that is, if and only if K st = 0, then the equality holds. Besides, it is assumed that the static equilibrium position is at the maximum negative stiffness of TMO as commented earlier, then we can conclude that TMOriVI has the minimum stiffness at the equilibrium position. Thus, the change trend of stiffness is opposite on both sides of the equilibrium position, namely, where > 0 describes the neighborhood of u st . Therefore, the static equilibrium position is an inflection point of the restoring force curve for TMOriVI. Then, we have On the other hand, combing the Taylor expansions Eq. (51), we can calculate the derivatives as follows In summary, we can obtain the following conclude: for the dynamic system described in Eq. (43) and Eq. (51), the Taylor series satisfies k 3 = 0, k 2 = 0, k 1 ≥ 0, and k 1 = 0 if and only if K st = 0.