Theoretical Analysis and Numerical Simulation of Resonances and Stability of a Piecewise Linear-Nonlinear Vibration Isolation System

A methodology is presented to study the resonance and stability for a single-degree-of-freedom (SDOF) system with a piecewise linear-nonlinear stiffness term (i.e., one piece is linear and the other is weakly nonlinear). Firstly, the exact response of the linear governing equation is obtained, and a modified perturbation method is applied to finding the approximate solution of the weakly nonlinear equation. Then, the primary and 1/2 subharmonic resonances are obtained by imposing continuity conditions and periodicity conditions. Furthermore, Jacobian matrix is derived to investigate the stability of resonance responses. Finally, the results of theoretical study are compared with numerical results, and a good agreement is observed.


Introduction
Vibration isolation and shock absorbing have always been a hot research topic in engineering practice.A solid and liquid mixture (SALiM) vibration isolator was developed to isolate vibrations and shocks induced by heavy machines [1][2][3].The stiffness of the SALiM isolator is piecewise linear-nonlinear in a quasi-static test [3].That is, the isolator exhibits linear stiffness in some displacement region, beyond which the nonlinear stiffness is observed.Therefore, the elastic restoring force is continuous, but its first-order derivative at the turning point is discontinuous.In the past two decades, the majority of researches focused on the dynamics response, stability, bifurcation, and chaos of piecewise linear systems, such as piecewise bilinear systems or piecewise trilinear systems [4][5][6][7][8].However, piecewise linear-nonlinear or even nonlinearnonlinear systems have not received much attention, but many physical systems in fields of aerospace engineering, electric circuit, and so forth appear to be piecewise linearnonlinear systems [9][10][11][12].
For nonsmooth stiffness systems, most approaches finding their steady state responses could be sorted into three groups including harmonic balance method (HBM) and its modified form, increment harmonic balance method (IHBM) [9,10], classical approximate analytical methods like average method [4,[13][14][15], and the matching method [6,7,11,16].For instance, using HBM, Jin et al. investigated the response and stability of an unsymmetrical multiple-degreeof-freedom system with piecewise linear elastic elements [9], and for a more complex oscillator possessing a periodically time-varying and piecewise binonlinear restoring force function HBM is still applicable to its periodic motion [9,17].But like all harmonic balance techniques, the accuracy of the results obtained by HBM and IHBM depends on the number of harmonics included.In practice, engineers prefer some classic approximate approaches like average method in engineering practice.That is because such techniques could provide a closed form of frequency response function which is useful for dynamics design.Hu studied the primary resonance of a harmonically forced oscillator with a pair of symmetric setup elastic stops using average approach [15] and developed design criteria [8].For suspension mechanisms, Jazar et al. completed frequency response analysis and jump avoidance formula for a suspension system using average method [18] and Deshpande et al. presented a comprehensive optimal design solution for a similar piecewise nonlinear suspension [14].However, the use of average method is confined by the presumption of weak nonlinearity, but piecewise systems are considered as strongly nonlinear due to the frequent occurrence of bifurcations that result in subharmonics or chaotic behaviors.
Both approaches mentioned above can be applied to approximate solutions, but more substantial and sophisticated dynamic properties induced by nonsmooth nonlinearity cannot be explored thoroughly.Matching method provides the exact form of periodic solution, but some unknown coefficients have to be determined by the numerical method.The return map technique associating with matching method is preferable to explore more detailed nonlinear phenomena involving stability analysis, bifurcations, and chaos [19][20][21], because the exact solution from matching method could give an adequately accurate Jacobian matrix, the eigenvalues of which determine the stability of steady state solutions.For instance, Shaw and Holmes gave a classic academic article and applied return map method to study a periodically forced piecewise linear oscillator on aspects of time-domain response, subharmonics, and chaos motions [7].Andreaus et al. modeled a mechanical impact system with soft contact using ordinary differential equation with piecewise bilinear terms and implemented a numerical iterative method in association with return map, based on which bifurcation diagrams and the Lyapunov coefficients were obtained [22].For a symmetric nonlinear-linear single-degree-of-freedom magnetic bearing system, Ji and Hansen studied symmetric periodic solution using the matching method [11].Chatterjee et al. transferred a harmonically excited Duffing oscillator with a one-sided elastic restraint to a piecewise linear system and gave the corresponding approximate periodic solution [12].
The present investigation intends to extend the previous studies by developing an analytical procedure which determines the single-crossing periodic responses of a class of harmonically excited piecewise linear-nonlinear oscillators.And also, the developed procedure can determine the crossing time as well as the stability of located periodic orbits.
The organization of this paper proceeds as follows.Firstly, the mathematical model is abstracted from a real bellows type SALiM isolator.And primary and subharmonic resonances are evaluated with matching method in association with perturbation method which is utilized to seek approximate solution of weakly nonlinear equation.Based on the preceding discussion, the Jacobian matrix is established with return map method and thus the stability of located periodic resonance orbits is determined, respectively, using the two different approaches.Finally, one or two examples will be introduced to demonstrate the validity of the approaches.

Equation of Motion
For the bellows type SALiM isolator, as the isolator is stretched at a quasi-static speed, there will come a point where elastic restoring force of the solid elements filled in the container is completely balanced by its internal pressure, which can be referred to as discontinuity point in a sense of stiffness [3].After that point, the volume of the elements will no longer increase as the isolator is further stretched.This means that at discontinuity point the stiffness of the isolator suddenly becomes softer, which is equal to the container's linear stiffness.Hence, the displacement corresponding to the stiffness discontinuity point can be referred to as critical displacement.The vibration isolation system discussed herein could be simplified into a single-degree-of-freedom massspring system with a piecewise linear-nonlinear stiffness as shown in Figure 1.
The equation of motion of the system in Figure 1 can be written as where  designates the displacement of gravity center of the mass ,  is damping coefficient,   () is elastic restoring force of the piecewise linear-nonlinear stiffness, and  and  are the amplitude and frequency of external excitation force.Equation (1) can be rewritten as where  1 = /,  1 = /, and the form of elastic restoring force is given as in [3] where  1 ,  2 , and  3 are stiffness coefficients of nonlinear piece and  is the linear stiffness coefficient.It should be noted that  0 is the coordinate value of discontinuity point and not an absolute distance.
Obviously, there exist two kinds of solutions for the piecewise linear-nonlinear system, namely, a small amplitude vibration (which means that the displacement does not exceed the critical displacement point) and a large amplitude vibration which can be beyond the critical point.The main work of this paper is to seek the large amplitude response of the overall system.

Resonance Analysis
3.1.Primary Resonance.This section presents a detailed analytical procedure to find approximate solutions for the overall system under primary and subharmonic resonances.The analysis is based on the assumption that there exists an asymptotic expansion solution for the nonlinear differential equation with a weak nonlinearity [11].In order to establish the solution for the large amplitude response of the overall system, it is necessary firstly to seek the individual general solutions of ordinary differential equations corresponding to linear and nonlinear stiffness, respectively.Then both solutions for the two different stiffness sections are matched at the discontinuity points of stiffness by applying the continuity conditions and periodicity conditions of displacement and velocity.According to the linear ordinary differential equation theory, the general solution including transient response induced by initial conditions and steady forced response is easy to obtain.However, the general solution for the nonlinear dynamic equation is not available so far, so an approximate solution is needed.Perturbation method is a commonly used approach and will be utilized below.
Before further discussion, (2) can be rewritten as two equations as The single-crossing steady state response is depicted in Figure 2, which shows that the overall response during one period consists of two segments in the following two time intervals: [ 0 ,  1 ] and [ 1 ,  2 ], where  0 ,  1 , and  2 denote the time instants at the critical displacement.The motion in each time interval is locally governed by the corresponding equations (4) and (5).Especially, let  1 () be the system displacement response from (4), starting at time  0 and ending at time  1 , and  2 () indicates the displacement response during the period [ 1 ,  2 ], which is the general solution of (5).
For the first segment of the motion, an approximate solution for the primary resonance response can be determined using an asymptotic expansion method.In interval [ 0 ,  1 ], the corresponding displacement response represented by the solution of (4) could be solved by considering the displacement and velocity at  0 as initial conditions.Without loss of generality, initial time  0 is assumed to be zero, but a new initial excitation phase  0 must be given as reparation.In order to employ a perturbation method, the following variables are introduced:  1 = ,  1 = ,  2 =  2 , and  3 =  3 .To avoid the small divisor terms in the first-order approximate solution in the vicinity of primary resonance, a detuning parameter is introduced according to  1 =  2 + , where  is the external detuning parameter.
It is assumed that an approximate periodic solution to (4) is a straightforward expansion of asymptotic series with respect to powers of the small parameter [11].The first-order approximate solution for ( 4) is given by  1 =  0 +  1 and thus ẋ 1 = u 0 +  u 1 and ẍ 1 = ü 0 +  ü 1 .Here the second and higher order terms are neglected in the asymptotic solution, although any desired higher order terms can be easily obtained using a similar procedure.
Substituting the above relating parameters and asymptotic solution into (4) and equating the coefficients of  0 and  1 yield It can be seen in ( 7) that excitation force, damping ratio, and two nonlinear terms appear.That is an indication that the first-order approximate solution can provide an adequate approximation of the exact solution, because the additional high order terms in the asymptotic solution make little contribution because of the small perturbation parameter.The solution of ( 6) is a well-known standard linear response: where  1 and  1 are constants to be determined.Then substituting ( 8) into (7) gives where the coefficients   ( = 1, 2, . . ., 7) will not be detailed for brevity.It is simple to obtain the solution of (9): Therefore, the solution to (4) could be directly obtained by combining ( 8) and (10).Actually, considering that two secular terms do not have enough time to grow and lead to an unbounded response in the steady state response, two secular terms need not be eliminated [11].That is because the time interval [ 0 ,  1 ] of this motion is limited and small.This point is different from conventional multiscale perturbation method.
The displacement response  2 () in time interval [ 1 ,  2 ] is considered below.According to the theory of linear ordinary differential equation, the solution of (5) with initial conditions at instant  1 is obtained: where  2 and  2 are constants to be determined and Till now, both of displacement histories  1 and  2 are obtained, but there are six unknown constants including  1 ,  1 ,  2 ,  2 , crossing time  1 , and initial excitation phase  0 in (8), (10), and (11).These parameters can be determined by imposing initial conditions, continuity and periodicity conditions: Substituting (19) and ( 20) into ( 13)-( 18) yields six equations with respect to  1 ,  1 ,  2 ,  2 ,  1 , and  0 .After a series of complicated algebra operations, two transcendental equations with respect to  1 and  0 are obtained.Therefore, the complete response history of the system cannot be obtained in closed form, and numerical means have to be adopted to find the crossing time  1 and initial phase  0 .
3.2.1/2 Subharmonic Resonance.The analytical procedures for seeking an approximate solution of 1/2 subharmonic resonance are implemented similarly.The transform of  =  +  2 /3 3 is introduced herein to further simplify the following process.Thus ( 4) and ( 5) can be rewritten as where In order to study 1/2 subharmonic resonance of the examined system, the following variables are introduced:  1 =  and  =  0 .Just as the discussion in Section 3.1, to avoid the appearance of the small divisor terms in the firstorder approximate solution in the vicinity of 1/2 subharmonic resonance, a detuning parameter is introduced according to  2  = (1/4) 2 + , where  is also the external detuning parameter.The approximate solution to the first order for ( 21) is  1 =  0 +  1 and thus ẏ 1 = u 0 +  u 1 and ÿ 1 = ü 0 +  ü 1 and the second and higher order terms are neglected in the asymptotic solution.
Abstracting the coefficients of  0 and  1 gives The solution of ( 23) is derived as where Ã1 and B1 are unknown constants for the homogeneous equation of (23) and Substituting ( 25) into (24) yields a complicated equation and again the coefficients   ( = 1, 2, . . ., 23) are not detailed for brevity: Although (26) appears to be complicated, it is simple to obtain its exact analytical solution according to the linear ordinary differential equation theory: Next, the solution of ( 24) can be obtained straightforwardly: where Again, the six constants, Ã , B ( = 1, 2),  1 , and  0 , need to be determined by imposing the initial conditions, continuity and periodicity conditions.The key point is that, for the 1/2 subharmonic resonance, the period must be altered as After obtaining   ( = 1,2), the inverse transform of  =  +  2 /3 3 gives the expected 1/2 subharmonic response,  1 and  2 .It should be emphasized that the methodology developed in this section is merely suitable for the single-crossing periodic response, which satisfies the form of solution depicted in Figure 2. The response for the multicrossing situation can also be obtained by extending the above procedure, at the expense of more cumbersome mathematics processing.

Return Map and Stability Analysis
In this section, the return map based on the Poincare section is employed to study the stability of periodic response.Generally, in order to establish the Poincare section for a singledegree-of-freedom periodically forced smooth oscillator, one frequently extends the two-dimensional phase plane (, ẋ ) to a three-dimensional phase space with coordinate variables (, ẋ , ) where time  is considered as the third state variable.The exact period of forced vibration of a nonautonomous system should be 2/, where  is the external excitation frequency.Thus for a smooth vibrating system, one often defines the Poincare section as [23][24][25] And this Poincare section is suitable for the nonsmooth system with piecewise linear-nonlinear stiffness as well.But the problem is that the original (undisturbed) periodic orbit and the nearby disturbed orbit starting from that section simultaneously cannot reach the discontinuous point where  =  0 at the same time.And hence a switching matrix (or saltation matrix) is introduced to sort out this problem in the process of establishing the Jacobian matrix [20,26].In fact, there exists an alternative Poincare section for the nonsmooth system [20,22], which is defined by Clearly, the fixed point in the Poincare section describes state variables ( ẋ , ) rather than (, ẋ ) which are utilized for the Poincare section (30).
4.1.Analytical Procedure.The three-dimensional space and the Poincare section represented by (31) are depicted in Figure 3.As can be seen, the starting point (denoted by ) with state variables (V 0 ,  0 ) in Σ  is mapped to another point (denoted by ) specified by (V 2 ,  2 ) lying in Σ  by the solution flow during one period.This process can be described by the map [7] P : Σ  → Σ  .
(32) And as can be seen in Figure 3, the entire space is divided into two subspaces by the section Σ  .Therefore, the map (32) is composed of two submaps: one is the map related to (4) and the other is governed by (5).The motion from point  to point  (specified by (V 1 ,  1 )) is governed by ( 4), and it is noted that point  is not in Σ  but in {(, ẋ , ) |  =  0 , ẋ < 0}.Therefore, the submap from points  to  is described by Similarly, the other submap from points  to  can be defined as Therefore, the complete map P can be expressed by a composition of the above submaps: with the Jacobian matrix where P 1 and P 2 denote the corresponding Jacobian matrixes of submaps P 1 and P 2 .
Next, the detailed procedures for seeking the Jacobian matrixes of both submaps will be studied.Firstly, considering that the state parameters (V 1 ,  1 ) are dependent on initial conditions (V 0 ,  0 ) and thus can be regarded as the functions of initial conditions, the Jacobian matrix of submap P 1 can be determined by At instant  0 , the initial conditions for the undisturbed trajectory are given as follows: where  1 () and ẋ 1 () have been obtained in Section 3.1.Note that the time-domain response  1 includes two unknown constants  1 and  1 , which are determined by initial conditions at initial instant  0 .Therefore, both  1 and  1 can be considered as the implicit functions of V 0 and  0 , and consequently initial conditions (38) could be rewritten as When the solution of (4) moves to point  at time  1 , the following continuity conditions hold: Differentiating (39) with respect to  0 and V 0 , respectively, gives Then differentiating (40) with respect to  0 and V 0 , respectively, yields Considering that the expressions of  1 and ẋ 1 are available in previous study,  1 / 1 ,  1 / 1 ,  ẋ 1 / 1 , and  ẋ 1 / 1 can be found directly.Further, one can obtain  1 /V 0 ,  1 / 0 ,  1 /V 0 , and  1 / 0 from the set of (43)-(44).
Similarly, taking derivatives of (41) with respect to  0 and V 0 gives Applying the same procedure to (42) gives Then after substituting  1 /V 0 ,  1 / 0 ,  1 /V 0 , and  1 / 0 obtained above into (45)-( 46), V 1 /V 0 , V 1 / 0 ,  1 /V 0 , and  1 / 0 could be obtained from linear equations (45)-(46).At this stage, the Jacobian matrix of submap P 1 can be constructed by the above four terms.Similarly, the Jacobian matrix of submap P 2 is expressed by In fact, P 2 can be obtained simply by replacing variables' subscripts "0" and "1" with "1" and "2" accordingly in the equation of P 1 , except the "0" in  0 .Once the two Jacobian matrices P 1 and P 2 are available, one can obtain P of the complete map (V 0 ,  0 ) → (V 2 ,  2 ) according to expression (36).Let  be the eigenvalue of matrix P with the larger modulus; if || < 1, it is indicated that the primary resonance response is stable.On the other hand, the solution is unstable when || > 1.

Numerical Scheme.
The previous analytical procedure is implemented based on the Poincare section defined by (31), which is frequently utilized to investigate how many times the periodic trajectory crosses the Poincare map.In order to compute numerically the Jacobian matrix, the Poincare section (30) is chosen.
In fact, there are two intersection points x 1 and x 2 of the border and the periodic orbit of primary resonance, as shown in Figure 4. Without losing generality, the neighborhood indicated by Σ 1 of x 1 is appointed as the Poincare section.And thus the corresponding map is described as Obviously, this section overlays that defined by (31), but the descriptions of state variables in the sections are distinguished.The whole map is composed by four submaps: where P 1 : Σ  are accordingly governed by the solution flows of ( 49) and (50) in subspaces  − and  + , respectively.And P 1 and P 3 are dominated by ( 49) and (50), respectively.P 2 and P 4 switch the state variables from one subspace to the other.Then how to seek the Jacobian matrixes of the four submaps will be explained.In subspace  − , considering that submap P 1 is smooth, the orbit starting from initial conditions x 1 in the Poincare section satisfies Taking the derivative of expression (53) with respect to x at x = x 1 yields Then the derivative of (54) with respect to time  gives Considering P 1 (x) =   (x) [27], the Jacobian matrix of submap P 1 satisfies the equation P 1 (x 1 ) =   10 (x 1 ), where subscript  10 denotes the time during which the orbit moves from Σ 1− to Σ 2− in subspace  − .Therefore, P 1 (x 1 ) = Φ( 10 ) where Φ() satisfying the initial condition of Φ( 0 ) = Φ(0) = I is the fundamental solution matrix of (49).Considering that (49) is nonlinear with respect to state variables, it is difficult to obtain Φ( 10 ) analytically.An alternative approach is a numerical method like the Runge-Kutta scheme to find the solution of the initial-value problem with x =  − (x, ) , Similarly, the Jacobian matrix P 2 (x 2 ) of the other submap P 2 may be calculated in the same way.Having established the Jacobian matrices of submaps P 1 and P 2 , the next step should be to find the Jacobian matrixes of P 2 and P 4 .Actually, P 2 is a local map representing the switch transform of state variables before and after the orbit crossing the border.And P 2 is not differentiable anywhere due to the discontinuity at the break point.In fact, P 2 (named as saltation matrix or switch matrix) symbolizes the relationship of variation x + to x − , which denote the variations between the undisturbed orbit x() and disturbed x() near the border, as can be seen in Figure 5 [27,28].In retrospect, P 1 and P 3 essentially demonstrate the variations of disturbed and undisturbed orbits, respectively, in subspaces  − and  + .The cause that generates the saltation matrix can be explained by the fact that the original and disturbed orbits cannot reach the border simultaneously as is depicted in Figure 5.At time   , original orbit x(  ) reaches the border Σ to a turn, but the disturbed orbit x(  ) lags behind by a short distance to the border.

Shock and Vibration
When the orbit goes across the boundary from  − to  + , the switching matrix is [28] where ∇ℎ is the gradient direction of the border defined in (48).Exchanging the subscript minus for a plus gives P 4 .Especially, the original orbit and disturbed orbit will return to the section Σ 1 at the same time and then set forth again, because both orbits possess the same periodicity and thus P 4 actually is an identity matrix and can be omitted.Finally,

Comparisons between Analytical and Numerical Results
To validate the approximate approaches, the periodic solutions determined by the developed analysis method are compared with the solutions from the numerical integration method.The classical fourth-order fixed-step Runge-Kutta algorithm is employed for the numerical integration.It is found that the approximate solutions and the numerical solutions are in excellent agreement for the cases of primary resonance and 1/2 subharmonic resonance.
As an example, the system's parameters are given as follows:  = 952,  1 = 1120,  2 = −3.53×10 4 ,  3 = 7.03×10      1 = 4.3,  1 = 0.5, and  0 = −0.001.Figure 6 shows the maximum displacement frequency response.The approximate solution is in good agreement with the numerical reference.In detail, in the case of excitation frequency  = 30 rad/s, the approximate analysis method in Section 3.1 finds  0 = 4.955939,  1 = 0.129850,  1 = 0.002196,  2 = −0.001894, 1 = −0.001178,and  2 = −0.008232.Figure 7 plots the periodic orbits of analytical approximate solution and solution obtained with numerical method.The maximum relative error between analytical displacement and numerical reference is 7.37% and is considered to be acceptable.Because the solution of nonlinear equation is just the result from approximate analysis in Section 3, it could explain the difference between periodic orbits in Figure 7.In this plot, the dot  line indicates the approximate analysis solution related to linear stiffness and the circle line represents the response of the weakly nonlinear equation.For the same system, the stability analysis developed in Section 4.1 gives a pair of eigenvalues with modulus |0.6382 ± 0.2945| = 0.70, which proves that the periodic orbit is stable.This conclusion could be confirmed by the result of |0.6039 ± 0.2040| = 0.64 obtained from numerical methodology developed in Section 4.2.The analytical procedure of finding eigenvalues is based on the approximate periodic response, and therefore the slight difference of eigenvalues from both methods might be explained principally by that reason.
To examine the 1/2 subharmonic resonance, consider a system with parameters  = 400,  1 = 1120,  2 = −3.53×10 4 ,  3 = 7.03 × 10 5 ,  1 = 2.13,  1 = 4,  0 = −0.001,and  = 56.55rad/s.The procedure discussed in Section 3.2 gives  0 = 6.212398,  1 = 0.092883, Ã1 = 0.004015, Ã2 = −0.002287,B1 = 0.000885, and B2 = −0.007041.Figure 8 shows the corresponding phase trajectory from approximate analysis method.For comparison, the numerical result is shown in the same plot.Apparently, the approximate solution agrees very well with the numerical result and the accuracy of the approximations should be adequate.The discrepancies are caused by the first-order truncation of the solution expansion of nonlinear equation (4).FFT technique is implemented in responses to inspect the frequency components.The related frequency spectrums are depicted in Figure 9.It manifests that zero frequency, primary frequency (approximating to 1/2 excitation frequency), and excitation frequency appear and other frequency components may be ignored due to their little contribution.By comparison, the peak amplitudes of the analytical response are slightly larger than those from

Figure 2 :
Figure 2: Periodic large amplitude response of the overall system.

Figure 3 :
Figure 3: The three-dimensional trajectory and the Poincare section.

Figure 4 :
Figure 4: The two-dimensional periodic orbit and the Poincare map.

Figure 5 :
Figure 5: The original and disturbed orbits.

Figure 6 :
Figure 6: Comparison between the displacement amplitudes from approximation analysis and numerical method.

Figure 7 :
Figure 7: Comparison of the periodic orbits from approximate analysis and numerical methods at  = 30 rad/s.

Figure 8 :Figure 9 :
Figure 8: Comparison of the periodic orbits of 1/2 subharmonic resonance response between the approximate and numerical solutions at  = 56.55rad/s.