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

Authors: X. Gao; Q. Chen
Journal: Shock and Vibration
 (2014)

---

## Abstract

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.


---

## Body

## 1. 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–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–8]. However, piecewise linear-nonlinear or even nonlinear-nonlinear systems have not received much attention, but many physical systems in fields of aerospace engineering, electric circuit, and so forth appear to be piecewise linear-nonlinear systems [9–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–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-degree-of-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–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.


## 2. 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 mass-spring system with a piecewise linear-nonlinear stiffness as shown in Figure 1.Figure 1
A single-degree-of-freedom mass-spring system.The equation of motion of the system in Figure1 can be written as
(1)Mx¨+Cx˙+Fk(x)=Fcos⁡ωt,
where x designates the displacement of gravity center of the mass M, C is damping coefficient, Fk(x) is elastic restoring force of the piecewise linear-nonlinear stiffness, and F and ω are the amplitude and frequency of external excitation force.Equation (1) can be rewritten as
(2)x¨+ξ1x˙+fk(x)=F1cos⁡(ωt),
where ξ1=C/M, F1=F/M, and the form of elastic restoring force is given as in [3]
(3)fk(x)={α1x+α2x2+α3x3(x≥a0)β(x-a0)+α1a0+α2a02+α3a03(x<a0),
where α1, α2, and α3 are stiffness coefficients of nonlinear piece and β is the linear stiffness coefficient. It should be noted that a0 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.


## 3. 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
(4)x¨+ξ1x˙+α1x+α2x2+α3x3=F1cos⁡(ωt+ϕ0)hhhhhhhhhhhhhhhhhhhhhhhhhwithx≥a0,(5)x¨+ξ1x˙+k2(x-a0)+α1a0+α2a02+α3a03=F1cos⁡(ωt+ϕ0)withx<a0.The single-crossing steady state response is depicted in Figure2, which shows that the overall response during one period consists of two segments in the following two time intervals: [t0,t1] and [t1,t2], where t0, t1, and t2 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 x1(t) be the system displacement response from (4), starting at time t0 and ending at time t1, and x2(t) indicates the displacement response during the period [t1,t2], which is the general solution of (5).Figure 2
Periodic large amplitude response of the overall system.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[t0,t1], the corresponding displacement response represented by the solution of (4) could be solved by considering the displacement and velocity at t0 as initial conditions. Without loss of generality, initial time t0 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=εξ, F1=εf, α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 x1=u0+εu1 and thus x˙1=u˙0+εu˙1 and x¨1=u¨0+εu¨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
(6)u¨0+ω2u0=0,(7)u¨1+ω2u1=-ξu˙0-σu0-α-2u02-α-3u03+fcos⁡(ωt+ϕ0).
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:
(8)u0=A1sin(ωt)+B1cos⁡(ωt),
where A1 and B1 are constants to be determined. Then substituting (8) into (7) gives
(9)u¨1+ω2u1=P1sin(ωt)+P2cos⁡(ωt)+P3sin(2ωt)+P4cos⁡(2ωt)+P5sin(3ωt)+P6cos⁡(3ωt)+fcos⁡(ωt+ϕ0)+P7,
where the coefficients Pi (i=1,2,…,7) will not be detailed for brevity. It is simple to obtain the solution of (9):
(10)u1(t)=-12P1tcos⁡(ωt)ω+12P2tsin(ωt)ω-13P3sin(2ωt)ω2-13P4cos⁡(2ωt)ω2-18P5sin(3ωt)ω2-18P6cos⁡(3ωt)ω2+P7ω2+14f[2ωtsin(ωt+ϕ0)+cos⁡(ωt+ϕ0)]ω2.
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 [t0,t1] of this motion is limited and small. This point is different from conventional multiscale perturbation method.The displacement responsex2(t) in time interval [t1,t2] is considered below. According to the theory of linear ordinary differential equation, the solution of (5) with initial conditions at instant t1 is obtained:
(11)x2=e-(1/2)ξ1(t-t1)[A2cos⁡ωd(t-t1)+B2sinωd(t-t1)]+[h1sin(ωt+ϕ0)-h2cos⁡(ωt+ϕ0)],
where A2 and B2 are constants to be determined and
(12)ωd=12-ξ12+4β,h1=F1ωξ1ω4+(ξ12-2β)ω2+β2,h2=F1(-β+ω2)ω4+(ξ12-2β)ω2+β2.Till now, both of displacement historiesx1 and x2 are obtained, but there are six unknown constants including A1, B1, A2, B2, crossing time t1, and initial excitation phase ϕ0 in (8), (10), and (11). These parameters can be determined by imposing initial conditions, continuity and periodicity conditions:
(13)x1(t0)=a0,(14)x1(t1)=a0,(15)x2(t1)=a0,(16)x2(t2)=a0,(17)x˙1(t1)=x˙2(t1),(18)x˙1(t0)=x˙2(t2),(19)t0=0,(20)t2=t0+2πω=2πω.Substituting (19) and (20) into (13)–(18) yields six equations with respect to A1, B1, A2, B2, t1, and ϕ0. After a series of complicated algebra operations, two transcendental equations with respect to t1 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 t1 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 ofy=x+α2/3α3 is introduced herein to further simplify the following process. Thus (4) and (5) can be rewritten as
(21)y¨+ξ1y˙+ωn2y+αy3=F1cos⁡(ωt+ϕ0)+Fc1,x>a-0,(22)y¨+ξ1y˙+βy=F1cos⁡(ωt+ϕ0)+Fc2,x<a-0,
where ωn2=α1-α22/3α3, α=α3, Fc1=α1α2/3α3-2α23/27α32, and Fc2=α2β/3α3-[(α1-β)a0+α2a02+α3a03]. And hence the critical displacement is also transferred into a-0=a0+α2/3α3.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 first-order approximate solution in the vicinity of 1/2 subharmonic resonance, a detuning parameter is introduced according to ωn2=(1/4)ω2+εσ, where σ is also the external detuning parameter. The approximate solution to the first order for (21) is y1=u0+εu1 and thus y˙1=u˙0+εu˙1 and y¨1=u¨0+εu¨1 and the second and higher order terms are neglected in the asymptotic solution.Abstracting the coefficients ofε0 and ε1 gives
(23)u¨0+14ω2u0=F1cos⁡(ωt+ϕ0)+Fc1,(24)u¨1+14ω2u1=-ξ1u˙0-σu0-α0u03.
The solution of (23) is derived as
(25)u0=A~1sin(ω2t)+B~1cos⁡(ω2t)+C,
where A~1 and B~1 are unknown constants for the homogeneous equation of (23) and C=-(4F1/3ω2)cos⁡(ωt+ϕ0)+4Fc1/ω2.Substituting (25) into (24) yields a complicated equation and again the coefficients Pi (i=1,2,…,23) are not detailed for brevity:
(26)u¨1+14ω2u1=P1sin(12ωt)+P2cos⁡(12ωt)+P3sin(32ωt)+P4cos⁡(32ωt)+P5sin(12ωt+ϕ0)+P6cos⁡(12ωt+ϕ0)+P7sin(32ωt+ϕ0)+P8cos⁡(32ωt+ϕ0)+P9sin(32ωt+2ϕ0)+P10cos⁡(32ωt+2ϕ0)+P11sin(2ωt+ϕ0)+P12cos⁡(2ωt+ϕ0)+P13cos⁡⁡(2ωt+2ϕ0)+P14cos⁡⁡(3ωt+3ϕ0)+P15sin⁡(52ωt+2ϕ0)+P16cos⁡(52ωt+2ϕ0)+P17sin(ωt)+P18cos⁡(ωt)+P19sin(ϕ0)+P20cos⁡(ϕ0)+P21sin(ωt+ϕ0)+P22cos⁡(ωt+ϕ0)+P23.Although (26) appears to be complicated, it is simple to obtain its exact analytical solution according to the linear ordinary differential equation theory:
(27)u1=-P1tωcos⁡(ωt2)+P2tωsin(ωt2)-P32ω2sin(3ωt2)-P42ω2cos⁡(3ωt2)+P5[sin(ωt/2+ϕ0)-ωtcos⁡(ωt/2+ϕ0)]ω2+P6[cos⁡(ωt/2+ϕ0)+ωtsin(ωt/2+ϕ0)]ω2-P72ω2sin(3ωt2+ϕ0)-P82ω2cos⁡(3ωt2+ϕ0)-P92ω2sin(3ωt2+2ϕ0)-P102ω2cos⁡(3ωt2+2ϕ0)-415P11sin⁡(2ωt+ϕ0)ω2-415P12cos⁡⁡(2ωt+ϕ0)ω2-415P13cos⁡(2ωt+2ϕ0)ω2-435P14cos⁡(3ωt+3ϕ0)ω2-16P15sin((5/2)ωt+2ϕ0)ω2-16P16cos⁡((5/2)ωt+2ϕ0)ω2-43P17sin(ωt)ω2-43P18cos⁡(ωt)ω2+4P19sin(ϕ0)ω2+4P20cos⁡(ϕ0)ω2-43P21sin(ωt+ϕ0)ω2-43P22cos⁡⁡(ωt+ϕ0)ω2+4P23ω2.Next, the solution of (24) can be obtained straightforwardly:
(28)y2=e-(1/2)ξ1(t-t1)[A~2cos⁡ωd(t-t1)+B~2sinωd(t-t1)]+[h3sin(ωt+ϕ0)-h4cos⁡(ωt+ϕ0)]+Fc2k2,
where h3=F1ωξ1/(ω4+(ξ12-2β)ω2+k22), h4=F1(-β+ω2)/(ω4+(ξ12-2β)ω2+β2), and ωd=(1/2)-ξ12+4β.Again, the six constants,A~i, B~i (i=1,2), t1, 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
(29)t2=t0+4πω=4πω.After obtainingyi (i=1,2), the inverse transform of y=x+α2/3α3 gives the expected 1/2 subharmonic response, x1 and x2. 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.


## 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
(4)x¨+ξ1x˙+α1x+α2x2+α3x3=F1cos⁡(ωt+ϕ0)hhhhhhhhhhhhhhhhhhhhhhhhhwithx≥a0,(5)x¨+ξ1x˙+k2(x-a0)+α1a0+α2a02+α3a03=F1cos⁡(ωt+ϕ0)withx<a0.The single-crossing steady state response is depicted in Figure2, which shows that the overall response during one period consists of two segments in the following two time intervals: [t0,t1] and [t1,t2], where t0, t1, and t2 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 x1(t) be the system displacement response from (4), starting at time t0 and ending at time t1, and x2(t) indicates the displacement response during the period [t1,t2], which is the general solution of (5).Figure 2
Periodic large amplitude response of the overall system.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[t0,t1], the corresponding displacement response represented by the solution of (4) could be solved by considering the displacement and velocity at t0 as initial conditions. Without loss of generality, initial time t0 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=εξ, F1=εf, α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 x1=u0+εu1 and thus x˙1=u˙0+εu˙1 and x¨1=u¨0+εu¨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
(6)u¨0+ω2u0=0,(7)u¨1+ω2u1=-ξu˙0-σu0-α-2u02-α-3u03+fcos⁡(ωt+ϕ0).
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:
(8)u0=A1sin(ωt)+B1cos⁡(ωt),
where A1 and B1 are constants to be determined. Then substituting (8) into (7) gives
(9)u¨1+ω2u1=P1sin(ωt)+P2cos⁡(ωt)+P3sin(2ωt)+P4cos⁡(2ωt)+P5sin(3ωt)+P6cos⁡(3ωt)+fcos⁡(ωt+ϕ0)+P7,
where the coefficients Pi (i=1,2,…,7) will not be detailed for brevity. It is simple to obtain the solution of (9):
(10)u1(t)=-12P1tcos⁡(ωt)ω+12P2tsin(ωt)ω-13P3sin(2ωt)ω2-13P4cos⁡(2ωt)ω2-18P5sin(3ωt)ω2-18P6cos⁡(3ωt)ω2+P7ω2+14f[2ωtsin(ωt+ϕ0)+cos⁡(ωt+ϕ0)]ω2.
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 [t0,t1] of this motion is limited and small. This point is different from conventional multiscale perturbation method.The displacement responsex2(t) in time interval [t1,t2] is considered below. According to the theory of linear ordinary differential equation, the solution of (5) with initial conditions at instant t1 is obtained:
(11)x2=e-(1/2)ξ1(t-t1)[A2cos⁡ωd(t-t1)+B2sinωd(t-t1)]+[h1sin(ωt+ϕ0)-h2cos⁡(ωt+ϕ0)],
where A2 and B2 are constants to be determined and
(12)ωd=12-ξ12+4β,h1=F1ωξ1ω4+(ξ12-2β)ω2+β2,h2=F1(-β+ω2)ω4+(ξ12-2β)ω2+β2.Till now, both of displacement historiesx1 and x2 are obtained, but there are six unknown constants including A1, B1, A2, B2, crossing time t1, and initial excitation phase ϕ0 in (8), (10), and (11). These parameters can be determined by imposing initial conditions, continuity and periodicity conditions:
(13)x1(t0)=a0,(14)x1(t1)=a0,(15)x2(t1)=a0,(16)x2(t2)=a0,(17)x˙1(t1)=x˙2(t1),(18)x˙1(t0)=x˙2(t2),(19)t0=0,(20)t2=t0+2πω=2πω.Substituting (19) and (20) into (13)–(18) yields six equations with respect to A1, B1, A2, B2, t1, and ϕ0. After a series of complicated algebra operations, two transcendental equations with respect to t1 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 t1 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 ofy=x+α2/3α3 is introduced herein to further simplify the following process. Thus (4) and (5) can be rewritten as
(21)y¨+ξ1y˙+ωn2y+αy3=F1cos⁡(ωt+ϕ0)+Fc1,x>a-0,(22)y¨+ξ1y˙+βy=F1cos⁡(ωt+ϕ0)+Fc2,x<a-0,
where ωn2=α1-α22/3α3, α=α3, Fc1=α1α2/3α3-2α23/27α32, and Fc2=α2β/3α3-[(α1-β)a0+α2a02+α3a03]. And hence the critical displacement is also transferred into a-0=a0+α2/3α3.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 first-order approximate solution in the vicinity of 1/2 subharmonic resonance, a detuning parameter is introduced according to ωn2=(1/4)ω2+εσ, where σ is also the external detuning parameter. The approximate solution to the first order for (21) is y1=u0+εu1 and thus y˙1=u˙0+εu˙1 and y¨1=u¨0+εu¨1 and the second and higher order terms are neglected in the asymptotic solution.Abstracting the coefficients ofε0 and ε1 gives
(23)u¨0+14ω2u0=F1cos⁡(ωt+ϕ0)+Fc1,(24)u¨1+14ω2u1=-ξ1u˙0-σu0-α0u03.
The solution of (23) is derived as
(25)u0=A~1sin(ω2t)+B~1cos⁡(ω2t)+C,
where A~1 and B~1 are unknown constants for the homogeneous equation of (23) and C=-(4F1/3ω2)cos⁡(ωt+ϕ0)+4Fc1/ω2.Substituting (25) into (24) yields a complicated equation and again the coefficients Pi (i=1,2,…,23) are not detailed for brevity:
(26)u¨1+14ω2u1=P1sin(12ωt)+P2cos⁡(12ωt)+P3sin(32ωt)+P4cos⁡(32ωt)+P5sin(12ωt+ϕ0)+P6cos⁡(12ωt+ϕ0)+P7sin(32ωt+ϕ0)+P8cos⁡(32ωt+ϕ0)+P9sin(32ωt+2ϕ0)+P10cos⁡(32ωt+2ϕ0)+P11sin(2ωt+ϕ0)+P12cos⁡(2ωt+ϕ0)+P13cos⁡⁡(2ωt+2ϕ0)+P14cos⁡⁡(3ωt+3ϕ0)+P15sin⁡(52ωt+2ϕ0)+P16cos⁡(52ωt+2ϕ0)+P17sin(ωt)+P18cos⁡(ωt)+P19sin(ϕ0)+P20cos⁡(ϕ0)+P21sin(ωt+ϕ0)+P22cos⁡(ωt+ϕ0)+P23.Although (26) appears to be complicated, it is simple to obtain its exact analytical solution according to the linear ordinary differential equation theory:
(27)u1=-P1tωcos⁡(ωt2)+P2tωsin(ωt2)-P32ω2sin(3ωt2)-P42ω2cos⁡(3ωt2)+P5[sin(ωt/2+ϕ0)-ωtcos⁡(ωt/2+ϕ0)]ω2+P6[cos⁡(ωt/2+ϕ0)+ωtsin(ωt/2+ϕ0)]ω2-P72ω2sin(3ωt2+ϕ0)-P82ω2cos⁡(3ωt2+ϕ0)-P92ω2sin(3ωt2+2ϕ0)-P102ω2cos⁡(3ωt2+2ϕ0)-415P11sin⁡(2ωt+ϕ0)ω2-415P12cos⁡⁡(2ωt+ϕ0)ω2-415P13cos⁡(2ωt+2ϕ0)ω2-435P14cos⁡(3ωt+3ϕ0)ω2-16P15sin((5/2)ωt+2ϕ0)ω2-16P16cos⁡((5/2)ωt+2ϕ0)ω2-43P17sin(ωt)ω2-43P18cos⁡(ωt)ω2+4P19sin(ϕ0)ω2+4P20cos⁡(ϕ0)ω2-43P21sin(ωt+ϕ0)ω2-43P22cos⁡⁡(ωt+ϕ0)ω2+4P23ω2.Next, the solution of (24) can be obtained straightforwardly:
(28)y2=e-(1/2)ξ1(t-t1)[A~2cos⁡ωd(t-t1)+B~2sinωd(t-t1)]+[h3sin(ωt+ϕ0)-h4cos⁡(ωt+ϕ0)]+Fc2k2,
where h3=F1ωξ1/(ω4+(ξ12-2β)ω2+k22), h4=F1(-β+ω2)/(ω4+(ξ12-2β)ω2+β2), and ωd=(1/2)-ξ12+4β.Again, the six constants,A~i, B~i (i=1,2), t1, 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
(29)t2=t0+4πω=4πω.After obtainingyi (i=1,2), the inverse transform of y=x+α2/3α3 gives the expected 1/2 subharmonic response, x1 and x2. 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.


## 4. 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 single-degree-of-freedom periodically forced smooth oscillator, one frequently extends the two-dimensional phase plane(x,x˙) to a three-dimensional phase space with coordinate variables (x,x˙,t) where time t 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–25]
(30)Σt={(x,x˙,t)∣t=mod⁡(2πω)}.
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 x=a0 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
(31)Σx={(x,x˙,t)∣x=a0,x˙>0}.
Clearly, the fixed point in the Poincare section describes state variables (x˙,t) rather than (x,x˙) 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 A) with state variables (v0,t0) in Σx is mapped to another point (denoted by C) specified by (v2,t2) lying in Σx by the solution flow during one period. This process can be described by the map [7]
(32)P:Σx⟶Σx.Figure 3
The three-dimensional trajectory and the Poincare section.And as can be seen in Figure3, the entire space is divided into two subspaces by the section Σx. 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 A to point B (specified by (v1,t1)) is governed by (4), and it is noted that point A is not in Σx but in {(x,x˙,t)∣x=a0,x˙<0}. Therefore, the submap from points A to B is described by
(33)P1:(v0,t0)⟶(v1,t1).
Similarly, the other submap from points B to C can be defined as
(34)P2:(v1,t1)⟶(v2,t2).
Therefore, the complete map P can be expressed by a composition of the above submaps:
(35)P=P2∘P1
with the Jacobian matrix
(36)DP=DP2∘DP1,
where DP1 and DP2 denote the corresponding Jacobian matrixes of submaps P1 and P2.Next, the detailed procedures for seeking the Jacobian matrixes of both submaps will be studied. Firstly, considering that the state parameters (v1,t1) are dependent on initial conditions (v0,t0) and thus can be regarded as the functions of initial conditions, the Jacobian matrix of submap P1 can be determined by
(37)DP1=∂(v1(v0,t0),t1(v0,t0))∂(v0,t0)=[∂v1∂v0∂v1∂t0∂t1∂v0∂t1∂t0].
At instant t0, the initial conditions for the undisturbed trajectory are given as follows:
(38)a0=x1(t0),v0=x˙1(t0),
where x1(t) and x˙1(t) have been obtained in Section 3.1. Note that the time-domain response x1 includes two unknown constants A1 and B1, which are determined by initial conditions at initial instant t0. Therefore, both A1 and B1 can be considered as the implicit functions of v0 and t0, and consequently initial conditions (38) could be rewritten as
(39)a0=x1(A1(v0,t0),B1(v0,t0),t0),(40)v0=x˙1(A1(v0,t0),B1(v0,t0),t0).
When the solution of (4) moves to point B at time t1, the following continuity conditions hold:
(41)a0=x1(v0,t0,t1(v0,t0),A1(v0,t0),B1(v0,t0)),(42)v1(v0,t0)=x˙1(v0,t0,t1(v0,t0),A1(v0,t0),B1(v0,t0)).
Differentiating (39) with respect to t0 and v0, respectively, gives
(43)∂a0∂t0=∂x1∂A1∂A1∂t0+∂x1∂B1∂B1∂t0+∂x1∂t0=0,∂a0∂v0=∂x1∂A1∂A1∂v0+∂x1∂B1∂B1∂v0=0.
Then differentiating (40) with respect to t0 and v0, respectively, yields
(44)∂a0∂t0=∂x˙1∂A1∂A1∂t0+∂x˙1∂B1∂B1∂t0+∂x˙1∂t0=0,∂v0∂v0=∂x˙1∂A1∂A1∂v0+∂x˙1∂B1∂B1∂v0=1.
Considering that the expressions of x1 and x˙1 are available in previous study,∂x1/∂A1, ∂x1/∂B1, ∂x˙1/∂A1, and ∂x˙1/∂B1 can be found directly. Further, one can obtain ∂A1/∂v0, ∂A1/∂t0, ∂B1/∂v0, and ∂B1/∂t0 from the set of (43)-(44).Similarly, taking derivatives of (41) with respect to t0 and v0 gives
(45)∂a0∂t0=∂x1∂A1∂A1∂t0+∂x1∂B1∂B1∂t0+∂x˙1∂t1∂t1∂t0+∂x˙1∂t0=0,∂a0∂v0=∂x1∂A1∂A1∂v0+∂x1∂B1∂B1∂v0+∂x1∂t1∂t1∂v0+∂x1∂v0=0.
Applying the same procedure to (42) gives
(46)∂v1∂t0=∂x˙1∂A1∂A1∂t0+∂x˙1∂B1∂B1∂t0+∂x˙1∂t1∂t1∂t0+∂x˙1∂t0,∂v1∂v0=∂x˙1∂A1∂A1∂v0+∂x˙1∂B1∂B1∂v0+∂x˙1∂t1∂t1∂v0+∂x˙1∂v0.
Then after substituting ∂A1/∂v0, ∂A1/∂t0, ∂B1/∂v0, and ∂B1/∂t0 obtained above into (45)-(46), ∂v1/∂v0, ∂v1/∂t0, ∂t1/∂v0, and ∂t1/∂t0 could be obtained from linear equations (45)-(46). At this stage, the Jacobian matrix of submap P1 can be constructed by the above four terms.Similarly, the Jacobian matrix of submapP2 is expressed by
(47)DP2=∂(v2(v1,t1),t2(v1,t1))∂(v1,t1)=[∂v2∂v1∂v2∂t1∂t2∂v1∂t2∂t1].In fact,DP2 can be obtained simply by replacing variables’ subscripts “0’’ and “1’’ with “1’’ and “2’’ accordingly in the equation of DP1, except the “0’’ in a0. Once the two Jacobian matrices DP1 and DP2 are available, one can obtain DP of the complete map (v0,t0)→(v2,t2) according to expression (36). Let λ be the eigenvalue of matrix DP 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.


### 4.2. 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.As described before, the whole space could be divided into two separated subspaces{S-∣h(x)<0} and {S+∣h(x)>0} by the border
(48){Σ∣h(x)=x-a0=0}.
Actually, the complete periodic orbit lies in four parts: {S-∣h(x)<0}, {S+∣h(x)>0}, {Σ∣h(x)=0,x˙>0}, and {Σ∣h(x)=0,x˙<0}. Actually, the vibrating system with piecewise linear-nonlinear stiffness can be described by the equations of the following forms:
(49)x=f-(t,x)ifx∈S-,(50)x=f+(t,x)ifx∈S+,
where f- and f+ denote state equations, respectively, with nonlinear and linear stiffness and x is the state vector including displacement and velocity.In fact, there are two intersection pointsxb1 and xb2 of the border and the periodic orbit of primary resonance, as shown in Figure 4. Without losing generality, the neighborhood indicated by Σ1 of xb1 is appointed as the Poincare section. And thus the corresponding map is described as
(51)P:Σ1⟶Σ1.Figure 4
The two-dimensional periodic orbit and the Poincare map.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:
(52)P=P4∘P3∘P2∘P1,
where P1:Σ1-→Σ2-, P2:Σ2-→Σ2+, P3:Σ2+→Σ1+, and P4:Σ1+→Σ1-. The points in Σi- and Σi+ (i=1,2) are accordingly governed by the solution flows of (49) and (50) in subspaces S- and S+, respectively. And P1 and P3 are dominated by (49) and (50), respectively. P2 and P4 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 subspaceS-, considering that submap P1 is smooth, the orbit starting from initial conditions xb1 in the Poincare section satisfies
(53)φt(x)=x+∫0tf-(φs(x))ds.
Taking the derivative of expression (53) with respect to x at x=xb1 yields
(54)Dφt(xb1)=I+∫0tDf-(φs(xb1))Dφs(xb1)ds.
Then the derivative of (54) with respect to time t gives
(55)ddtDφt(xb1)=Df-(φt(xb1))Dφt(xb1)=A(t)Dφt(xb1).
Considering P1(x)=φt(x) [27], the Jacobian matrix of submap P1 satisfies the equation DP1(xb1)=Dφt10(xb1), where subscript t10 denotes the time during which the orbit moves from Σ1- to Σ2- in subspace S-. Therefore, DP1(xb1)=Φ(t10) where Φ(t) satisfying the initial condition of Φ(t0)=Φ(0)=I is the fundamental solution matrix of (49). Considering that (49) is nonlinear with respect to state variables, it is difficult to obtain Φ(t10) analytically. An alternative approach is a numerical method like the Runge-Kutta scheme to find the solution of the initial-value problem
(56)Φ˙=Df-(x,t)Φ,Φ(0)=I
with
(57)x=f-(x,t),x(0)=xb1.
Similarly, the Jacobian matrix DP2(xb2) of the other submap P2 may be calculated in the same way.Having established the Jacobian matrices of submapsP1 and P2, the next step should be to find the Jacobian matrixes of P2 and P4. Actually, P2 is a local map representing the switch transform of state variables before and after the orbit crossing the border. And P2 is not differentiable anywhere due to the discontinuity at the break point. In fact, DP2 (named as saltation matrix or switch matrix) symbolizes the relationship of variation δxb+ to δxb-, which denote the variations between the undisturbed orbit x(t)and disturbed x~(t) near the border, as can be seen in Figure 5 [27, 28]. In retrospect, DP1 and DP3 essentially demonstrate the variations of disturbed and undisturbed orbits, respectively, in subspaces S- and S+. 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 tb, original orbit x(tb) reaches the border Σ to a turn, but the disturbed orbit x~(tb) lags behind by a short distance to the border.Figure 5
The original and disturbed orbits.When the orbit goes across the boundary fromS- to S+, the switching matrix is [28]
(58)DP2=I+(f+-f-)·∇h∇h·f-,
where ∇h is the gradient direction of the border defined in (48). Exchanging the subscript minus for a plus gives DP4. 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 DP4 actually is an identity matrix and can be omitted. Finally,
(59)DP=DP4∘DP3∘DP2∘DP1=DP3∘DP2∘DP1.


## 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 A) with state variables (v0,t0) in Σx is mapped to another point (denoted by C) specified by (v2,t2) lying in Σx by the solution flow during one period. This process can be described by the map [7]
(32)P:Σx⟶Σx.Figure 3
The three-dimensional trajectory and the Poincare section.And as can be seen in Figure3, the entire space is divided into two subspaces by the section Σx. 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 A to point B (specified by (v1,t1)) is governed by (4), and it is noted that point A is not in Σx but in {(x,x˙,t)∣x=a0,x˙<0}. Therefore, the submap from points A to B is described by
(33)P1:(v0,t0)⟶(v1,t1).
Similarly, the other submap from points B to C can be defined as
(34)P2:(v1,t1)⟶(v2,t2).
Therefore, the complete map P can be expressed by a composition of the above submaps:
(35)P=P2∘P1
with the Jacobian matrix
(36)DP=DP2∘DP1,
where DP1 and DP2 denote the corresponding Jacobian matrixes of submaps P1 and P2.Next, the detailed procedures for seeking the Jacobian matrixes of both submaps will be studied. Firstly, considering that the state parameters (v1,t1) are dependent on initial conditions (v0,t0) and thus can be regarded as the functions of initial conditions, the Jacobian matrix of submap P1 can be determined by
(37)DP1=∂(v1(v0,t0),t1(v0,t0))∂(v0,t0)=[∂v1∂v0∂v1∂t0∂t1∂v0∂t1∂t0].
At instant t0, the initial conditions for the undisturbed trajectory are given as follows:
(38)a0=x1(t0),v0=x˙1(t0),
where x1(t) and x˙1(t) have been obtained in Section 3.1. Note that the time-domain response x1 includes two unknown constants A1 and B1, which are determined by initial conditions at initial instant t0. Therefore, both A1 and B1 can be considered as the implicit functions of v0 and t0, and consequently initial conditions (38) could be rewritten as
(39)a0=x1(A1(v0,t0),B1(v0,t0),t0),(40)v0=x˙1(A1(v0,t0),B1(v0,t0),t0).
When the solution of (4) moves to point B at time t1, the following continuity conditions hold:
(41)a0=x1(v0,t0,t1(v0,t0),A1(v0,t0),B1(v0,t0)),(42)v1(v0,t0)=x˙1(v0,t0,t1(v0,t0),A1(v0,t0),B1(v0,t0)).
Differentiating (39) with respect to t0 and v0, respectively, gives
(43)∂a0∂t0=∂x1∂A1∂A1∂t0+∂x1∂B1∂B1∂t0+∂x1∂t0=0,∂a0∂v0=∂x1∂A1∂A1∂v0+∂x1∂B1∂B1∂v0=0.
Then differentiating (40) with respect to t0 and v0, respectively, yields
(44)∂a0∂t0=∂x˙1∂A1∂A1∂t0+∂x˙1∂B1∂B1∂t0+∂x˙1∂t0=0,∂v0∂v0=∂x˙1∂A1∂A1∂v0+∂x˙1∂B1∂B1∂v0=1.
Considering that the expressions of x1 and x˙1 are available in previous study,∂x1/∂A1, ∂x1/∂B1, ∂x˙1/∂A1, and ∂x˙1/∂B1 can be found directly. Further, one can obtain ∂A1/∂v0, ∂A1/∂t0, ∂B1/∂v0, and ∂B1/∂t0 from the set of (43)-(44).Similarly, taking derivatives of (41) with respect to t0 and v0 gives
(45)∂a0∂t0=∂x1∂A1∂A1∂t0+∂x1∂B1∂B1∂t0+∂x˙1∂t1∂t1∂t0+∂x˙1∂t0=0,∂a0∂v0=∂x1∂A1∂A1∂v0+∂x1∂B1∂B1∂v0+∂x1∂t1∂t1∂v0+∂x1∂v0=0.
Applying the same procedure to (42) gives
(46)∂v1∂t0=∂x˙1∂A1∂A1∂t0+∂x˙1∂B1∂B1∂t0+∂x˙1∂t1∂t1∂t0+∂x˙1∂t0,∂v1∂v0=∂x˙1∂A1∂A1∂v0+∂x˙1∂B1∂B1∂v0+∂x˙1∂t1∂t1∂v0+∂x˙1∂v0.
Then after substituting ∂A1/∂v0, ∂A1/∂t0, ∂B1/∂v0, and ∂B1/∂t0 obtained above into (45)-(46), ∂v1/∂v0, ∂v1/∂t0, ∂t1/∂v0, and ∂t1/∂t0 could be obtained from linear equations (45)-(46). At this stage, the Jacobian matrix of submap P1 can be constructed by the above four terms.Similarly, the Jacobian matrix of submapP2 is expressed by
(47)DP2=∂(v2(v1,t1),t2(v1,t1))∂(v1,t1)=[∂v2∂v1∂v2∂t1∂t2∂v1∂t2∂t1].In fact,DP2 can be obtained simply by replacing variables’ subscripts “0’’ and “1’’ with “1’’ and “2’’ accordingly in the equation of DP1, except the “0’’ in a0. Once the two Jacobian matrices DP1 and DP2 are available, one can obtain DP of the complete map (v0,t0)→(v2,t2) according to expression (36). Let λ be the eigenvalue of matrix DP 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.


## 4.2. 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.As described before, the whole space could be divided into two separated subspaces{S-∣h(x)<0} and {S+∣h(x)>0} by the border
(48){Σ∣h(x)=x-a0=0}.
Actually, the complete periodic orbit lies in four parts: {S-∣h(x)<0}, {S+∣h(x)>0}, {Σ∣h(x)=0,x˙>0}, and {Σ∣h(x)=0,x˙<0}. Actually, the vibrating system with piecewise linear-nonlinear stiffness can be described by the equations of the following forms:
(49)x=f-(t,x)ifx∈S-,(50)x=f+(t,x)ifx∈S+,
where f- and f+ denote state equations, respectively, with nonlinear and linear stiffness and x is the state vector including displacement and velocity.In fact, there are two intersection pointsxb1 and xb2 of the border and the periodic orbit of primary resonance, as shown in Figure 4. Without losing generality, the neighborhood indicated by Σ1 of xb1 is appointed as the Poincare section. And thus the corresponding map is described as
(51)P:Σ1⟶Σ1.Figure 4
The two-dimensional periodic orbit and the Poincare map.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:
(52)P=P4∘P3∘P2∘P1,
where P1:Σ1-→Σ2-, P2:Σ2-→Σ2+, P3:Σ2+→Σ1+, and P4:Σ1+→Σ1-. The points in Σi- and Σi+ (i=1,2) are accordingly governed by the solution flows of (49) and (50) in subspaces S- and S+, respectively. And P1 and P3 are dominated by (49) and (50), respectively. P2 and P4 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 subspaceS-, considering that submap P1 is smooth, the orbit starting from initial conditions xb1 in the Poincare section satisfies
(53)φt(x)=x+∫0tf-(φs(x))ds.
Taking the derivative of expression (53) with respect to x at x=xb1 yields
(54)Dφt(xb1)=I+∫0tDf-(φs(xb1))Dφs(xb1)ds.
Then the derivative of (54) with respect to time t gives
(55)ddtDφt(xb1)=Df-(φt(xb1))Dφt(xb1)=A(t)Dφt(xb1).
Considering P1(x)=φt(x) [27], the Jacobian matrix of submap P1 satisfies the equation DP1(xb1)=Dφt10(xb1), where subscript t10 denotes the time during which the orbit moves from Σ1- to Σ2- in subspace S-. Therefore, DP1(xb1)=Φ(t10) where Φ(t) satisfying the initial condition of Φ(t0)=Φ(0)=I is the fundamental solution matrix of (49). Considering that (49) is nonlinear with respect to state variables, it is difficult to obtain Φ(t10) analytically. An alternative approach is a numerical method like the Runge-Kutta scheme to find the solution of the initial-value problem
(56)Φ˙=Df-(x,t)Φ,Φ(0)=I
with
(57)x=f-(x,t),x(0)=xb1.
Similarly, the Jacobian matrix DP2(xb2) of the other submap P2 may be calculated in the same way.Having established the Jacobian matrices of submapsP1 and P2, the next step should be to find the Jacobian matrixes of P2 and P4. Actually, P2 is a local map representing the switch transform of state variables before and after the orbit crossing the border. And P2 is not differentiable anywhere due to the discontinuity at the break point. In fact, DP2 (named as saltation matrix or switch matrix) symbolizes the relationship of variation δxb+ to δxb-, which denote the variations between the undisturbed orbit x(t)and disturbed x~(t) near the border, as can be seen in Figure 5 [27, 28]. In retrospect, DP1 and DP3 essentially demonstrate the variations of disturbed and undisturbed orbits, respectively, in subspaces S- and S+. 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 tb, original orbit x(tb) reaches the border Σ to a turn, but the disturbed orbit x~(tb) lags behind by a short distance to the border.Figure 5
The original and disturbed orbits.When the orbit goes across the boundary fromS- to S+, the switching matrix is [28]
(58)DP2=I+(f+-f-)·∇h∇h·f-,
where ∇h is the gradient direction of the border defined in (48). Exchanging the subscript minus for a plus gives DP4. 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 DP4 actually is an identity matrix and can be omitted. Finally,
(59)DP=DP4∘DP3∘DP2∘DP1=DP3∘DP2∘DP1.


## 5. 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×104, α3=7.03×105, ξ1=4.3, F1=0.5, and a0=-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, t1=0.129850, A1=0.002196, A2=-0.001894, B1=-0.001178, and B2=-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.2945i|=0.70, which proves that the periodic orbit is stable. This conclusion could be confirmed by the result of |0.6039±0.2040i|=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.Figure 6
Comparison between the displacement amplitudes from approximation analysis and numerical method.Figure 7
Comparison of the periodic orbits from approximate analysis and numerical methods atω=30 rad/s.To examine the 1/2 subharmonic resonance, consider a system with parametersβ=400, α1=1120, α2=-3.53×104, α3=7.03×105, ξ1=2.13, F1=4, a0=-0.001, and ω=56.55rad/s.The procedure discussed in Section 3.2 gives ϕ0=6.212398, t1=0.092883, A~1=0.004015, A~2=-0.002287, B~1=0.000885, and B~2=-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 numerical method and the peak frequencies coincide with each other perfectly.Figure 8
Comparison of the periodic orbits of 1/2 subharmonic resonance response between the approximate and numerical solutions atω=56.55 rad/s.Figure 9
Comparison between frequency spectrums from approximate analysis and numerical integration (○: numerical, □: analytical).Actually, the 1/2 subharmonic motion can also be observed by the bifurcation diagram. The period-doubling bifurcation indicates the existence of 1/2 subharmonic motion. The fourth-order Runge-Kutta numerical integration is performed to plot the bifurcation diagram, and the Poincare section used here is defined by (30). Taking the excitation frequency as the bifurcation parameter gives the bifurcation diagrams of the velocity response as shown in Figure 10. It is shown that the period-doubling motion appears in the frequency band between 8.28 and 9.98 Hz due to the loss of the stability of period-one motion and then vanishes again at 9.98 Hz.The bifurcation diagrams.
(a)
(b)The appearance or disappearance of period-doubling motion may be explained by the change of eigenvalues of the Jacobian matrix of period-one motion. For the inverse period-doubling bifurcation near 9.98 Hz, Figure11 gives the locations of eigenvalues of 9.98 and 10.0 Hz in the unit circle. It can be seen that all eigenvalues locate in the real axis and the most important fact is that the eigenvalue of larger modulus passes the unit circle by −1, which indicates the occurrence of inverse period-doubling bifurcation. In other words, the period-doubling motion turns into a period-one orbit when the excitation frequency shifts from 9.98 to 10.0 Hz. It is worth mentioning that the fact that the larger eigenvalue of 9.98 Hz lies outside of unit circle (see Figure 11) indicates that the period-one orbit is unstable, but period-doubling motion is stable, which has been exhibited in the bifurcation diagram of velocity in Figure 10. As further explanation, Figure 12 presents the phase orbits at 9.98 Hz and 10.00 Hz. The zoomed curves illustrate that a period-doubling motion exists at the frequency of 9.98 Hz, but no evidence supports the existence of the period-doubling motion at 10.0 Hz.Figure 11
Unit circle and eigenvalues of period-one and period-doubling motions. The eigenvalue with larger magnitude passes by −1 into unit circle from −1.022 to −0.9647.Phase orbits at different excitation frequencies: (a) 9.98 Hz and (b) 10.00 Hz.
(a)
(b)


## 6. Conclusions

In the first part of this paper, the matching method in association with the perturbation technique has been proven an efficient approach to find primary resonance and 1/2 subharmonic resonance responses. And in the second part, the return map based on the Poincare sections defined byΣt={(x,x˙,t)∣t=mod⁡(2π/ω)} and Σx={(x,x˙,t)∣x=a0,x˙>0} was employed to study the stability of periodic responses. The former section was applied to establish the Jacobian matrix of periodic resonance response analytically, and the numerical scheme was developed based on the latter definition. The consistence between the analytical and numerical eigenvalues of the periodic orbit verifies the stability analysis. Moreover, the developed numerical method can also be conveniently extended to the bifurcation analysis of periodic orbits. From the perspective of vibration isolation, the occurrence of any subharmonic motion is a threat to the performance of a vibration isolator and thus is supposed to be eliminated in the design stage.



---

*Source: 803275-2014-03-02.xml*