An archetypal zero- or quasi-zero-stiffness model with three degrees of freedom based upon an inverse method

In this paper, a three degrees of freedom archetypal model with arbitrary designable features and its inverse construction method based upon the target restoring force functions are proposed. The construction conditions of the zero-stiffness and quasi-zero-stiffness systems with 3-DOF have been studied in detail, by introducing the definitions of the zero-stiffness and the quasi-zero-stiffness system, based upon the archetypal model. The static characteristic and the low-frequency vibration isolation performance of the zero-stiffness and quasi-zero-stiffness systems with 3-DOF in multiple direction have been analyzed, by comparing with the linear system with the same support stiffness. The conclusion shows that there are finite solutions for the construction of 3-DOF zero-stiffness characteristic, and infinite solutions for the construction of 3-DOF quasi-zero-stiffness characteristic, based on 3-DOF archetypal model and its dimensionless parameters. The feasibility of constructing the zero-stiffness system with 3-DOF will lead to a new inspiration for the design of the passive zero-stiffness support structures; and the low-frequency vibration isolation performance of the 3-DOF quasi-zero-stiffness system in multiple direction is significantly improved compared with the linear system, which expresses a great potential in the multi-DOF design of the low-frequency, ultra-low-frequency vibration isolation.


Introduction
In many engineering fields represented by building aseismic [1][2][3][4][5], low-frequency vibration isolation [6][7][8], and the simulation of the free-boundary condition [9][10][11][12] or the ground simulation of zero or microgravity environment [13][14][15][16], there is an extremely strong demand for multi-degree-of-freedom (multi-DOF) support system with high-static-low-dynamic stiffness and even zero-natural-frequency (zero-stiffness). Since the quasi-zero-stiffness (QZS) system is proposed by Molyneux [17] in 1957 firstly, the QZS system has solved the paradox between the low-frequency vibration isolation and high bearing capacity in the field of low-frequency isolator design, based upon its unique characteristics of the high-static-low-dynamic stiffness. The complex dynamics behavior of the QZS system has been gradually revealed and deeply studied by large number of scholars [18][19][20][21][22][23]. At the same time, it has also been widely used in many low-frequency vibration isolation fields and made remarkable achievements [24][25][26]. Ibrahim [27] applied the principle of positive and negative stiffness in the low-frequency vibration isolation fields to seismograph firstly in 1974. Vaneijk et al. [28] constructed a leaf spring with QZS characteristics based on the negative stiffness charac-teristics of Euler compression bar structure at the University of Colombo in 1976. Carrella et al. [30] studied the application in the field of vibration isolation to reduce the whirling of a rotating shaft by simplifying the QZS model into a Duffing system with thirdorder approximation in 2009. Xu et al. [29] designed a QZS vibration isolator based upon magnet connecting rod at Hunan University in 2013. Wang et al. [31] designed an experimental prototype to study a compact QZS isolator using wave springs in 2021. It is noteworthy that some scholars have realized the design of vibration isolator with quasi-zero stiffness by using cam-roller mechanism, and proved its advantageous low-frequency vibration isolation performance through theoretical analysis, simulation calculation, and experimental researches. Figure 1 shows two examples of the quasi-zero stiffness vibration isolator based upon camroller mechanism, in which Fig. 1a is the QZS vibration isolator with cam-roller-spring mechanism designed by Zhou et al. [32], and Fig. 1b is the novel integrated quasi-zero stiffness vibration isolator for coupled translational and rotational vibrations designed by Ye et al. [33].
Although the nonlinear vibration isolation system with QZS shows great application value in the field of low-frequency vibration isolation, most of the research still focuses on single-DOF or can only show significant low-frequency vibration isolation performance in single-DOF. Therefore, in order to further broaden the application direction of QZS system in the field of lowfrequency vibration isolation engineering, many scholars have made exploratory research on the application of QZS system in the field of multi-DOF vibration isolation. Wang et al. [34] developed an active vibration isolation platform with 6-DOF based upon eight pneumatic actuators with spring-damping system in each cavity, to improve the vibration isolation performance even when the actuator fails. Sun et al. [35] designed a low-frequency vibration isolator with three translational DOF with QZS in the vertical direction by paralleling four horizontally arranged hinge mechanisms with vertical springs in 2015. Zhu et al. [36] designed a 6-DOF vibration isolation device with QZS in the vertical direction by using the magnetic levitation system combined with the control algorithm. Shi et al. [37] designed a low-frequency vibration isolator with QZS characteristics in the horizontal plane by paralleling two groups of QZS unit in 2016. Zhou et al. [38] designed a vibration isolation platform with 6-DOF and QZS characteristics based on Stewart mechanism by introducing 12 permanent magnet QZS compression rods in 2017. Xiao et al. [39] achieved 6-DOF vibration isolation performance by introducing a kind of QZS compression unit into Stewart platform. Zheng et al. [40] designed a parallel vibration isolation platform with 3-DOF QZS characteristics by using horizontal spring linkage instead of translational active pair in 2018. The simulation results show that the system has a significant improvement in low-frequency vibration isolation performance within 3-DOF. Zhu et al. [41] designed a QZS vibration isolator with 2-DOF in the horizontal directions based upon the two-stage vibration isolation method, to isolate the seismic wave for electrical equipments.
However, there are still some deficiencies in the design of multi-DOF nonlinear low-frequency vibration isolator. For example, some schemes inevitably need to introduce control systems to ensure the stability of the structure. And some schemes design the QZS system with multi-DOF by introducing the whole quasizero stiffness mechanism as an elastic element into the parallel vibration isolation system, which leads the system structure extremely complex, increasing the difficulty of the theoretical analysis and implementation greatly [42]. Therefore, the research of a multi-DOF passive low-frequency vibration isolation system based upon a relatively simple structure and strong applicability is still very necessary.
The motivation of this paper is to propose an 3-DOF archetypal model with arbitrary design characteristics and its inverse construction method based upon the target restoring force functions. Another motivation is to study the construction condition of the zero-stiffness (ZS) or the quasi-zero-stiffness (QZS) characteristic based upon the archetypal model, and then analyze the low-frequency vibration isolation performance of the QZS system with 3-DOF compared with the linear system with the same support stiffness.

Archetypal model with multiple degrees of freedom
The 3-DOF archetypal model shown in Fig. 2 comprises a mass block with cube shape, three groups of springs, three groups of irregular surfaces (which is highlighted in blue in Fig. 2), and three groups of guiding mechanisms with orthogonal characteristic (which Fig. 1 Examples of the quasi-zero stiffness vibration isolator based upon cam-roller mechanism, a for the QZS vibration isolator with cam-roller-spring mechanism designed by Zhou et al. [32], b for the novel integrated quasi-zero stiffness vibration isolator for coupled translational and rotational vibrations designed by Ye et al. [33] Gravity Direction  are highlight in red, green, and yellow in Fig. 2, respectively). Each group of springs consists of two linear springs with the same stiffness which are arranged collinear, respectively, and ensure that each spring can only be stretched or compressed along its axis without deflection. The axes of any two groups of springs are orthogonal to each other, and the intersection coincides with the center of gravity of the load under the condition of static equilibrium. The cosine of the angle between the axis of any group of springs and the gravity direction is √ 3/3. Each group of planes (a plane and its opposite) of the mass block is perpendicular to the axis of each group of springs, and a vertical axis is matched with the mass block along the gravity direction through a pair of linear bearings (as shown in the green part on the mass block in Fig. 2), to ensure that the relative displacement between the mass block and the vertical axis is always along the axis direction of the vertical axis. There is a pair of linear bearing (as shown in the yellow part on the vertical axis) assembled at both the ends of the vertical axis which is matched with a pair of horizontal axis, respectively, to ensure that the rel-ative displacement between the vertical axis and the horizontal axis is always along the axis direction of the horizontal axis. Both the ends of the each horizontal axis is matched with the longitudinal slide rail through the longitudinal slider (as shown in the red part on the horizontal axis), to ensure that relative displacement between the horizontal axis and the longitudinal slide rail is always along the axis direction of the longitudinal slide rail. Assuming the center of gravity of the mass block is located at the geometric center, and three groups of irregular surfaces are fixedly installed on the x-, y-, and z-axes planes of the block, and connected with the free ends of the six springs through six balls. Assuming each ball is always in frictionless contact with the surface within the effective range. The coordinate system X OY is established, by taking center of gravity of the load as the origin at the static equilibrium state; the X -axis is along the axis of one group of springs and the positive direction is specified as the inverse direction of the projection of the gravity vector on the axis; the Y -and Z -axis are determined according to the right-hand rule, and ensure that the projection of the gravity vector on the Y -or Z -axis is along its negative direction, respectively, as shown in Fig. 2.
Assuming the mass of the load is m; the total stiffness along X -, Y -, or Z -axes is k X , k Y , or k Z , respectively; the equilibrium length of the spring in positive or negative direction of the i-axis is L i+ or L i− , respectively, in which i = X , Y , or Z . Since the equilibrium length of the spring is designable, we stipulate that the equilibrium length of each the spring in X -, Y -, or Zaxis is always designed according to the following rela- in which g X = g Y = g Z = √ 3/3g ≈ 5.66m/s 2 express the component of the gravity in the negative direction of X -, Y -, and Z -axis, respectively; L i and H i (in which i = X , Y , or Z ) can be considered as the equilibrium length and the initial compression of each the spring equipped in i-axis under the condition of no gravity, respectively, in order to simplify the symbol. The three groups of irregular surfaces S X, SY , and S Z with everywhere differentiable feature are assumed to have rotational symmetry about the X -, Y -, and Z -axes, respectively, which can be defined as following in order to ensure the symmetry of the mechanical characteristics of the system. The origin of the coordinate system of each surface function is coincides with the center of gravity of the load, and the positive direction of the coordinate axis of the dependent variable of each surface is always in the same direction as the vector from the origin to the zero-point of the surface function and meets the right-hand rule.
The compression of each the spring can be expressed as following in which i+ or i− expresses the compression of the spring located at the positive or negative direction of the i-axis, and when the displacement of the load is (X, Y, Z ). Then, the elastic force vector generated by each the spring can be obtained as following The normal vector of each the ball contacts with the surface S X, SY , or S Z can be expressed as following in which Then, the reaction force vector provided by each spring to the load can be expressed as following Obviously, the restoring force vector of the system can be obtained as following in which Thus, the equation of motion can be expressed as following in which Bx 0 = Sx 0 − H X , By 0 = Sy 0 − H Y , and to the displacement R X , R Y , or R Z . Then, the dimensionless equation of motion can be written as the following , and s z (ρ z ) = Sz (R Z );ẋ,ẏ,ż,ẍ,ÿ, orz represents the first and the second derivatives of the dimensionless displacements x, y, or z to the dimensionless time τ , respectively. The surface functions of the multi-DOF archetypal model can be determined by solving the following differential equations where f x (x, y, z), f y (x, y, z), or f z (x, y, z) represents the component of the target restoring force function vector in the x-, y-, or z-axis direction, respectively.
After equivalent transformation (13) can be rewritten as following Then, the solution of the differential equations (13) can be obtained as following in which where C x , C y , or C z is the undetermined constant, which need to be determined by introducing the boundary condition of the target restoring force vector. Thus, the three groups of surfaces in the archetypal model with multi-DOF can always be solved uniquely and definitively by (15), based upon the target restoring force vector of the system and an appropriate boundary conditions.

Zero-stiffness condition with multi-DOF
This section will study the zero-stiffness condition with multi-DOF based upon the archetypal model with multi-DOF. Firstly, the definition of zero-stiffness system with multi-DOF is proposed as following where f (x) is a smooth function vector dependents on the dimensionless displacement vector x ∈ R n , which implies the equivalent stiffnessk (x) = df (x) dx T ≡ 0, |x| ∈ (0, 1) denoted as ZS n system. According to Definition 3.1, if condition (17) holds, the following equation must holds based upon the restoring force vector (14). Divide the sum of the three formulas in (18) by 2, and the following formula can be obtained Making the difference between (19) and each the formulas in (18), and the following equations can be obtained Since the differential equations (20) are decoupled, taking the first equation as an example consider the boundary value problem of the following differential equation when Rewrite s x (ρ x ) and s x (ρ x ) as s x and ds x /dρ x , respectively, and separate the variables s x and ρ x to both sides of the first formula in (21), then it can be rewritten as Performing indefinite integral on both sides of (22), there is Substituting the second formula of (21) into (23), and the undetermined constant C x can be determined as following Thus, the solution of the boundary value problem of the differential equation (21) can be obtained as following According to the physical meaning of the parameters κ x , κ y , and κ z , there is which is easy to prove α x = 0. It can be proved that when α x > 0, the surface s x = s x (ρ x ) defined by equation (25) is an ellipsoid, otherwise it is a hyperboloid with the focus on x-axis. It is easy to draw a conclusion by analogy: when α y > 0 or α z > 0, the surface s y = s y ρ y or s z = s z (ρ z ) defined by equation (25) is an ellipsoid, otherwise it is a hyperboloid with the focus on y-axis or z-axis. According to the mapping relationship between the parameters α x , α y , and α z and the parameters κ x , κ y , and κ z , it is easy to judge that there is at most one of the surface functions s x = s x (ρ x ), s y = s y ρ y , and s z = s z (ρ z ) may be hyperboloid. Therefore, the parameter space κ x , κ y , κ z can be divided into four regions, which are defined as following which represent the different combinations of three pairs of surface functions in the archetypal model with multi-DOF, and its distribution is shown in Fig. 3. When the parameters combination κ x , κ y , κ z in the red, green, blue, or yellow region, respectively, it indicates that the following four situations of the surfaces combination: 1. all the surfaces are ellipsoids (E 3 ); 2. the surface s x = s x (ρ x ) is hyperboloid and the surfaces s y = s y ρ y and s z = s z (ρ z ) are ellipsoids (H x ); 3. the surface s y = s y ρ y is hyperboloid and the surfaces s x = s x (ρ x ) and s z = s z (ρ z ) are ellipsoids (H y ); 4. the surface s z = s z (ρ z ) is hyperboloid and the surfaces s x = s x (ρ x ) and s y = s y ρ y are ellipsoids (H z ).
Next, consider the situation of the parameters combination κ x , κ y , κ z is not inside any color region but on the interface of any two regions, shown in Fig. 3. According to the geometric relationship, it is easy to prove that when the parameter combination falls on surface S x , S y , or S z , it means that the parameters separately meet the relationship as following or Taking (28) situation as an example to analyze the shapes of the surfaces, according to the symmetry characteristics of the system. The equations can be obtained as following in which β y = κ y /κ z and β z = κ z /κ y = 1/β y , by making a difference between (19) and each the formulas in (18). Firstly, consider the first formula in (31), which can derive the boundary value problem of the following differential equation as following According to the first formula in (32), there is Substituting the second formula in (32) into (33), the undetermined constant C x can be determined as following Substitute (34) into (33), and the following relationship can be obtained based upon the arbitrariness of the parameter s x0 . It is means that s x = s x (ρ x ) is a constant function, that is a plane perpendicular to the x-axis. Since forms of the second and the third formulas are completely consistent with the first formula of (21), the following conclusions can be drawn directly It is easy to know that β y > 0 and β z > 0 is always held, which means that both the surfaces identified by (36) are ellipsoids, based upon the definitions of β y and β z . In particular, if the stiffness in the y-axis and z-axis are equal, the two designable surfaces perpendicular to the y-axis and z-axis are spherical, and the radius should be equal to the initial compression h y and h z of the springs, respectively.

Dynamics analysis of the ZS 3 system
Based upon the uniqueness of the restoring force function of the zero-stiffness system, it is advisable to take a symmetric structure with three pairs of congruent ellipsoid surfaces as an example to analyze the dynamic characteristics of the ZS 3 system. Assuming in which where α > 0 and b are undetermined constants. Then, the equation of motion of system (12) can be rewritten as following The analytical transmissibilityfrequency relationships of the ZS 3 system, the quasizero stiffness vibration isolation system based on the Duffing system [43], and the linear system with the same supporting stiffness can be obtained easily as following and by introducing the ZS 3 condition into (39), and introducing the same periodic excitation function f (τ ) = f cos ωτ into these systems, respectively. The analytical solution of the transmissibility of the ZS, QZS, and linear systems which are plotted in red, blue, and black, respectively, in Fig. 4. The solid, dashed, and dotted-dashed represent the transmissibility with the damping ratio satisfy ζ = 0.1, 0.2, and 0.3, respectively. It is obviously that the QZS system can produce vibration attenuation effect at least above ω * QZS ≈ 0.56, while the ZS system can reduce the initial isolation frequency to ω * ZS = 0. It is widely known that the increase in the damping ratio will improve the low-frequency isolation capability of the system by suppressing the resonance phenomenon, at the expense of high-frequency isolation performance. And the demarcation point between the high-and lowfrequency section of a system can be considered as the initial vibration frequency, deservedly. Hence, the global isolation performance of ZS system will only increase with the decrease in the damping ratio ζ , for its characteristic of ω * ZS = 0 which indicates that the ZS system does not have a low-frequency band and resonance phenomenon. Further, if the damping ratio satisfied ζ = 0, the system will not have any ability of energy transmission, to realize the immune effect of vibration in full frequency band.

Quasi-zero-stiffness condition with multi-DOF
Firstly, the definition of QZS system with multi-DOF is proposed by analogy with the definition of single-DOF, in order to study the conditions of the QZS system with multi-DOF for system (12).
in which f (x) is a smooth function vector dependents on the dimensionless displacement vector x ∈ R n denoted as QZS n system. In particular, if the system is a quasi-zero-stiffness system with single-DOF, the lower right corner mark can be omitted and marked as QZS system. Furthermore, the system is called stablequasi-zero-stiffness with n-DOF (denoted as SQZS n ) system if and only if in which δ ∈ R n is a vector with arbitrarily small module, otherwise, the system is called unstable-quasizero-stiffness with n-DOF (denoted as UQZS n ).
Thus, the condition of QZS 3 for system (12) can be expressed as the following proposition.

Proposition 4.1 The necessary and sufficient conditions for the archetypal model with QZS 3 characteristics is
in which C is an arbitrary constant.
Proof The dimensionless restoring force function vector of system (12) can be expressed as following The dimensionless equivalent stiffness function vector can be obtained as followinĝ According to the definition of coordinate system of the archetypal model, it is obviously that the static equilibrium position x 0=0 , and the dimensionless restoring force vector satisfies Thus, the condition of the QZS 3 for the system can be expressed as following based upon Definition 4.1, in which ρ = (x, y, z) T expresses the dimensionless displacement vector. Introducing the following transformation to eliminate the influence from the direction of the vector ρ when it converges to 0, in which α and β are arbitrary angle variables independent of ρ, respectively; and it is easy to prove that ρ = |ρ| is the module of the dimensionless displacement vector ρ. According to the everywhere differentiable and the symmetry characteristics, it is obviously that Substituting (46), (49), and (50) into (48), the following equation must be held in which s x2 = d 2 s x (ρ x ) , and . Based upon the arbitrariness of the angles α and β, the following relationship must be held in order to ensure the establishment of equation (51), which is equivalent to the following ⎛ Then, equation (54) must have nonzero solutions, by considering the physical meaning of the parameters of κ x , κ y , and κ z , which is equivalent to hold the following equation And the solution of (55) which is also the necessary condition of the QZS 3 for system (12) can be obtained as following On the other hand, substituting (56) into (46), and leading the module of the dimensionless displacement vector ρ = (x, y, z) T approach to 0, it is easy to obtain the following relationship which means that the sufficiency of the QZS 3 condition for system (12) is proved. Thus, Proposition 4.1 is proved.

Modeling and static analysis of the QZS 3 system
It is easy to prove that the QZS 3 condition can be held if replace the combination of the mass block and all surfaces in the archetypal model to a sphere with radius l 0 . According to the characteristics of the sphere, the distance from the center of gravity of the load block to any point on the surface is always the same, so the restoring force of the system is decoupled from the rotational DOFs of the mass block. It is means that the system has ZS 3 characteristic within the rotational DOFs, and the guiding mechanisms shown in Fig. 2 can be omitted to simplified the model with QZS 3 characteristic in the translational DOFs, shown in Fig. 5. Substituting the spherical function into (12), and the motion equation of the system shown in Fig. 5 can be obtained as following It is easy to solve the dimensionless parameters based upon function (58). Substituting (60) into (45), and assuming the dimensionless stiffness κ x , κ y , and κ z satisfy the relationship as following then the following equations can be solved based upon the dimensionless relationship Substituting (62) into (56), then the equations can be obtained as following and the QZS 3 condition of system (59) can be solved, based upon the assumption of the damping ratio equal in all directions, as following The motion equation of the QZS 3 system can be obtained as following in which the restoring force can be expressed as following Considering the dimensionless restoring force function (65), its direction may not strict same as the vector of the dimensionless displacement when the direction of the displacement is not along any axis. Thus, the decomposition of the restoring force vector along its displacement direction (radial component) f ρ can be expressed as following Consider the following ratio to analyze the proximity between the restoring force amplitude of the restoring force and its radial component. It is easy to obtained the following relationship by Taylor series expansion where O (i) represents the same order infinitesimal of i, which means that the restoring force and its radial component are so close that it is reasonable to approximate the restoring force of the original QZS 3 system with its radial component. Then, the amplitude of the restoring force function of the approximate system can be obtained by introducing (49) into the radial component (67) as following in which λ 2 x = cos 2 α + sin 2 α sin 2 β, λ 2 y = cos 2 α + sin 2 α cos 2 β, λ 2 z = sin 2 α, σ x = γ x / γ x + γ y + γ z , , σ y = γ y / γ x + γ y + γ z , and , σ z = γ z / γ x +γ y +γ z . And the equivalent stiffness of the approximate system can be solved as followinĝ Firstly, suppose all the springs along every axes are equal, which means The motion equation of the symmetric QZS 3 system can be obtained as following by substituting (72) into (65) in which the restoring force can be expressed as following The restoring force curves along the x-, y-, and zaxes are express the same SQZS characteristics simultaneously, which are shown in Fig. 6.
Furthermore, the distribution of the dimensionless restoring force vector in the space (x, y, z) is shown in Fig. 7, in order to study the stability of the system in the whole displacement space. The homocentric spheres with different colors separately identified multiple surfaces with the same displacement from the initial equilibrium position, in which the closer the displacement is to the origin, the closer the color of the sphere is to red; and the arrows express the direction of the dimensionless restoring force, while the magnitudes of the vectors from large to small are identified by the various colors from blue to red. It is easy to find that the restoring force vector of the system at any position in space has the trend of pointing to the center of the ball; the closer the position is to the origin, the smaller the amplitude of the restoring force is; and when the position is on the same sphere, the amplitude of restoring force will gradually increase as its position approaches any coordinate axis.
The decomposition of the restoring force vector along its displacement direction (radial component) f ρ , and the decomposition along the tangent direction of the displacement direction (tangential component) f t can separately be expressed as following The exploded view of the restoring force vector is shown in Fig. 8a, in which the arrows in red or blue represent the radial component or the tangential component of the restoring force vector, and the magnitude of each the vector amplitude is indicated by the length of the arrow. It can be found that the tangential component is not zero and its amplitude is much smaller than the radial component when the displacement is not on any coordinate axis or any line of l 1 : x = y = z, l 2 : −x = y = z, l 3 : x = −y = z, or l 4 : x = y = −z, otherwise, the tangential component is zero. The distribution of the tangential component is shown separately in Fig. 8b, to analyze the distribution under the condition of eliminating the interference of the radial component. It is easy to find that the direction of the tangential component in any quadrant is always along the sphere and pointed from a coordinate plane to the line l 1 , l 2 , l 3 or l 4 , and the closer to the straight line, the smaller the amplitude of the vector will be.
Then, the dimensionless restoring force function of the approximate system can be obtained by introducing (72) into the radial component (70) as following The restoring force curves of the approximate and the original systems separately along the direction of axis (take the x-axis as an example), line in coordinate plane (take the direction of α = 2/3 β = π/2 as an example), or an arbitrary line in the displacement space (take the direction of α = 2/3 β = 3/4 as an example) are shown in Fig. 9, to analyze the deviation of the restoring force between the approximate and the original systems, based upon (74) and (77). The black and red curves in Fig. 9 express the restoring force of the original and approximate systems, respectively. It can be found that the approximate system and the original system have a high coincidence in the whole definition domain of displacement, which means that the approximate analysis of the original system using the approximate system can maintain a high accuracy. Thus, the conclusion of (69) is proved intuitively, that it is reasonable to use the approximate system composed of the radial component instead of the restoring force vector to approximate analyze the original system (73), based upon the premise that the amplitude of the tangential component is much smaller than that of radial component.
In the other hand, there are obvious differences in the equivalent stiffness of the system, when the displacement along different directions, which means that it is necessary to further analyze the distribution of the maximum and minimum values of the equivalent stiffness of the system. The equivalent stiffness of the system can be obtained based upon (77) as followinĝ When the displacement amplitude ρ = 0 is a constant, the stationary point (α s , β s ) of the equivalent stiffnesŝ k ρ with respect to the angle α and β can be expressed by solving the equations in which and It is obviously that the 6 stationary points in the set V A represent the directions along the positive and negative directions of each axis, respectively; and the 8 stationary points in the set V B represent the directions along the 4 lines of l 1 , l 2 , l 3 , and l 4 , respectively. According the sufficient conditions to determine the maximum or minimum of multivariate function, letting Substituting all the elements in the set of stagnation points (80) into (83), and the relationships separately adapt the sets V A and V B can be obtained as following and It is easy to prove that the following relationships always hold if ρ = 0 and which means that the equivalent stiffness is the maximum or minimum separately when the displacement direction is along any coordinate axis or the displacement direction is along any line in l 1 , l 2 , l 3 , and l 4 , on the premise of the same displacement amplitude.

Symmetric QZS 3 system
Introducing a periodic excitation along an arbitrary direction into the approximate system to analyze the vibration isolation performance of the symmetric QZS 3 system, and the equation of motion can be obtained as following Assuming the velocity and displacement solutions are periodic, which can be indicated as following where ψ (τ ) = ωτ + φ (τ ), which expressed by two variables of the amplitude a and phase φ, that can be assumed to be slowly varying functions of dimensionless time τ . Thus, the first and the second derivatives of ρ to τ can be obtained as following Substituting (88) and (89) into (90), the formulas describing the derivative of amplitude and phase can be obtained as following in which According the slowly varying characteristics of the parameters a and φ to τ , let P (a, φ) sin ψ and P (a, φ) cos ψ approximately equal to the average value of the integral over a period, that is Substituting (92) into (93), the average equation of the parameters a and φ can be solved as following in which in which EllipticK (m) or EllipticE (m) denotes the complete elliptic integral of the first-or the second-kind with the elliptic module m.
Let ȧ = 0 φ = 0 , that derive singular point (ρ s , φ s ) that corresponds to the periodic solution of steady state of system (88) which satisfies Then, the dimensionless amplitude-frequency relationship can be solved by eliminating the phase φ as following The analytical and the numerical solution of the amplitude-frequency curves with different parameters are shown in Fig. 10, in which Fig. 10a or b expresses the curve with different directions of the excitation or the different damping ratios, respectively. The dashed, solid, or dotted-dashed curve in Fig. 10a separately represents the analytical solution at the situation ), or the bisector between y-axis and the vector v that is the vector w = √ 3 + 3 /6, 1/ 6 + 2 √ 3, 1/ 6 + 2 √ 3 (i.e., α = arccos 1/ 6 + 2 √ 3 β = arccos 1/ 5 + 2 √ 3 ); the dashed, dotteddashed, or solid curve in Fig. 10b represents the analytical solution when the damping ratio satisfies ζ = 0.10, 0.15, or 0.20, respectively; the symbol +, ×, or • represents the numerical solution corresponding to the analytical solution plotted by the dashed, dotted-dashed, or solid curve, respectively, and the symbols in red or blue are separately express the numerical solutions solved by the forward-or the backward-sweep process. Comparing the analytical and numerical results, it can be found that the two kinds of results show a high degree of agreement, which proves the accuracy of the analytical conclusions. Particularly, the bifurcation occurs in the low-frequency band shown in Fig. 10a, due to the nonlinear characteristics of the system in multiple direction, which can be expressed that the amplitude response increases along the large-amplitude branch and then falls down to the small-amplitude branch at the downward jumping point, denoted by the red-arrow directions, as the excitation frequency increases from ω = 0 gradually. In the other case, when the excitation frequency decreases from ω = 1 gradually, the amplitude response increases along the small-amplitude branch and then jumps up to the largeamplitude branch at the upward jumping point, denoted by the blue-arrow directions. Combined with the stiffness analysis results of the system, it is easy to notice that when the equivalent stiffness in the excitation direction is small the frequency of the resonance peak will shift left, while the peak value of the system amplitude response will be amplified due to its soft equivalent stiffness. However, whether the jumping phenomenon caused by the nonlinear characteristics of the system or the increase in the resonance peak caused by the soft equivalent stiffness can be suppressed by enlarging the damping ratio appropriately, which can be concluded from Fig. 10b. The force transmissibility of system (88) can be defined as following in which F max is the maximum of the reaction force function F r (ρ) = 2ζρ in order to further analyze the vibration isolation performance. Because the reaction force function F r (ρ) has monotonicity in respect of the displacement ρ, the maximum value of the reaction force can be derived with an approximately non-conservative estimate as following It can be solved that the initial vibration isolation frequencies of the system are about 0.46, 0.44, and 0.42, when the directions of excitation are along the y-axis, vector w, and vector v, respectively. Compared with the initial vibration isolation frequency of the linear system ω * Linear = √ 2 ≈ 1.41, it can produce at least about 67.5% reduction, and the resonance peak will also produce at least about 35.9% suppression. Figure 12 shows the comparison between the initial vibration isolation frequencies of the linear system and the QZS 3 system when the excitation along different directions. The red and blue curves show infimum and the supremum of the initial vibration isolation frequency when the excitation along the directions of yaxis and vector v, respectively; the black curve shows the initial vibration isolation frequency of the linear system; and the region with the varying color shows all the possible situations of the initial vibration isolation frequency of QZS 3 system. It can be seen that the initial vibration isolation frequencies of QZS 3 system have solutions only when the amplitude of the external excitation is less than about 2 (otherwise the amplitude response a may beyond the definition domain of the system displacement), based upon the numerical analysis results, since there is no analytical solution of equation (102). From the upper bound of the initial vibration isolation frequency curve of the system shown as the red curve in Fig. 12 (the displacement direction is along the y-axis), it can be found that under the condition that the system has the worst low-frequency vibration isolation performance, the QZS 3 system can still ensure expresses the effect of reducing for the initial vibration isolation frequency compared with linear system when the excitation amplitude is less than about 1.97. In other words, QZS 3 system has lower initial vibration isolation frequency compared with linear system in almost all situations, and the difference of the low-frequency vibration isolation performance along different directions is tiny.

Asymmetric QZS 3 system
Similar to the analysis method of stiffness distribution in Sect. 4.2, it is easy to prove that the directions of maximum equivalent stiffness are independent of the stiffness of each spring, that is, the directions of maximum equivalent stiffness are always along the directions of the coordinate axes. Therefore, the initial vibration isolation frequency of the amplitude-frequency relationship or transmissibility of the asymmetric QZS 3 system along each coordinate axis can be considered as the upper bound in multiple direction, and the analysis of the vibration isolation performance of the system in multiple direction can be simplified to the analysis of the vibration isolation performance of the system along each coordinate axis. The vibration isolation performance of the system with different stiffness along each axis can be analyzed by introducing a periodic excitation along an arbitrary direction into (65), and the motion equation of this forced system can be expressed as following Similarly, the amplitude-frequency relationship of system (103) can be solved as following Z x (a)+ Z y (a) + Z z (a) + ω 2 2 + (2ζ ω) 2  EllipticK λ 2 y a 2 λ 2 y a 2 − 1 Fig. 13 Amplitude-frequency relationship with different stiffness along each direction of axis for the asymmetric QZS 3 system, a for the parameters satisfying δ x = 0.5 and δ y = 0.5, b for the parameters satisfying δ x = 1.0 and δ y = 0.5, c for the parameters satisfying δ x = 1.5 and δ y = 0.5, d for the parameters satisfying δ x = 0.5 and δ y = 1.0, e for the parameters satisfying δ x = 1.0 and δ y = 1.0, f for the parameters satisfying δ x = 1.5 and δ y = 1.0, g for the parameters satisfying δ x = 0.5 and δ y = 1.5, h for the parameters satisfying δ x = 1.0 and δ y = 1.5, and i for the parameters satisfying δ x = 1.5 and δ y = 1.5. (Color figure online) missibility of the nonlinear system are always much lower than the linear one. The upper bound of the initial vibration isolation frequency separately along the direction of the x-, y-, or z-axis under the condition of the changing stiffness and excitation is shown in Fig. 15. It can be found from Fig. 15a-f that when the condition δ x − δ y − 1 = 0, δ x − δ y + 1 = 0, or δ x + δ y − 1 = 0 holds, the upper bound of the initial vibration isolation frequency of the QZS 3 system along the x-, y-, or z-axis direction is 0, which means that it represents zero-stiffness character-istic along the axis direction. From Fig. 15d-f, when δ x increasing, it can be found that the upper bound of the initial vibration isolation frequency separately has the characteristic of decreasing, decreasing to 0 first then increasing, or increasing slightly, along the x-, y-, or z-axis direction, if δ y is constant and δ y > 1; if δ y is constant and δ y ≤ 1, the upper bound of the initial vibration isolation frequency separately has the characteristic of decreasing to 0 first then increasing, increasing, or decreasing to 0 first then increasing slightly, along the x-, y-, or z-axis direction, respectively; Simi-

Conclusion
In this paper, we have proposed a multi-DOF archetypal model with arbitrary designable features and its inverse construction method based upon the target restoring force functions. The construction conditions of the ZS 3 and the characteristics of the surfaces in the archetypal model are studied in detail, based upon the definition of the multi-DOF zero-stiffness system. The definition of multi-DOF quasi-zero-stiffness system is proposed based upon the definition of single-DOF quasi-zerostiffness system. The construction conditions of QZS 3 based upon the archetypal model with multi-DOF are obtained, and the static and the dynamics characteristics along multiple direction of the system are analyzed. The results show that the zero-stiffness characteristic along an arbitrary direction can be realized by introduce plane, sphere, ellipsoid, or hyperboloid into the archetypal model, based upon the dimensionless parameters of the system. The difference of the vibration isolation performance of the QZS 3 system along different directions are tiny, and the worst performance appears in the direction along the coordinate axis. Compared with the linear system with the same supporting stiffness, the QZS 3 system expresses the reduction ratio of the initial vibration isolation frequency at least about 67.5%, and the suppression effect of the resonance peak at least about 35.9%. It is means that this novel multi-DOF model expresses a great potential in the multi-DOF design of the low-frequency, ultra-low-frequency vibration isolation, and even the passive zero-stiffness support structures.