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.
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.
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 , , and the form of elastic restoring force is given as in [3] where , , and are stiffness coefficients of nonlinear piece and is the linear stiffness coefficient. It should be noted that 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
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: and , where , , and 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 be the system displacement response from (4), starting at time and ending at time , and indicates the displacement response during the period , 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 , the corresponding displacement response represented by the solution of (4) could be solved by considering the displacement and velocity at as initial conditions. Without loss of generality, initial time is assumed to be zero, but a new initial excitation phase must be given as reparation. In order to employ a perturbation method, the following variables are introduced: , , , and . 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 , 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 and thus and . 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 and 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 and are constants to be determined. Then substituting (8) into (7) gives where the coefficients () 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 of this motion is limited and small. This point is different from conventional multiscale perturbation method.
The displacement response in time interval is considered below. According to the theory of linear ordinary differential equation, the solution of (5) with initial conditions at instant is obtained: where and are constants to be determined and
Till now, both of displacement histories and are obtained, but there are six unknown constants including , , , , crossing time , and initial excitation phase 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 , , , , , and . After a series of complicated algebra operations, two transcendental equations with respect to and 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 and initial phase .
3.2. Subharmonic Resonance
The analytical procedures for seeking an approximate solution of 1/2 subharmonic resonance are implemented similarly. The transform of is introduced herein to further simplify the following process. Thus (4) and (5) can be rewritten as where , , , and . And hence the critical displacement is also transferred into .
In order to study 1/2 subharmonic resonance of the examined system, the following variables are introduced: and . 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 , where is also the external detuning parameter. The approximate solution to the first order for (21) is and thus and and the second and higher order terms are neglected in the asymptotic solution.
Abstracting the coefficients of and gives The solution of (23) is derived as where and are unknown constants for the homogeneous equation of (23) and
Substituting (25) into (24) yields a complicated equation and again the coefficients () 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 , , and .
Again, the six constants, , (), , and , 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 (), the inverse transform of gives the expected 1/2 subharmonic response, and . 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 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 , where is the external excitation frequency. Thus for a smooth vibrating system, one often defines the Poincare section as [23–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 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 () in is mapped to another point (denoted by ) specified by () lying in by the solution flow during one period. This process can be described by the map [7]
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 ()) is governed by (4), and it is noted that point is not in but in . Therefore, the submap from points to is described by Similarly, the other submap from points to can be defined as Therefore, the complete map can be expressed by a composition of the above submaps: with the Jacobian matrix where and denote the corresponding Jacobian matrixes of submaps and .
Next, the detailed procedures for seeking the Jacobian matrixes of both submaps will be studied. Firstly, considering that the state parameters () are dependent on initial conditions () and thus can be regarded as the functions of initial conditions, the Jacobian matrix of submap can be determined by At instant , the initial conditions for the undisturbed trajectory are given as follows: where and have been obtained in Section 3.1. Note that the time-domain response includes two unknown constants and , which are determined by initial conditions at initial instant . Therefore, both and can be considered as the implicit functions of and , and consequently initial conditions (38) could be rewritten as When the solution of (4) moves to point at time , the following continuity conditions hold: Differentiating (39) with respect to and , respectively, gives Then differentiating (40) with respect to and , respectively, yields Considering that the expressions of and are available in previous study,, , , and can be found directly. Further, one can obtain , , , and from the set of (43)-(44).
Similarly, taking derivatives of (41) with respect to and gives Applying the same procedure to (42) gives Then after substituting , , , and obtained above into (45)-(46), , , , and could be obtained from linear equations (45)-(46). At this stage, the Jacobian matrix of submap can be constructed by the above four terms.
Similarly, the Jacobian matrix of submap is expressed by
In fact, can be obtained simply by replacing variables’ subscripts “0’’ and “1’’ with “1’’ and “2’’ accordingly in the equation of , except the “0’’ in . Once the two Jacobian matrices and are available, one can obtain of the complete map according to expression (36). Let be the eigenvalue of matrix with the larger modulus; if , it is indicated that the primary resonance response is stable. On the other hand, the solution is unstable when .
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 and by the border Actually, the complete periodic orbit lies in four parts: , , , and . Actually, the vibrating system with piecewise linear-nonlinear stiffness can be described by the equations of the following forms: where and denote state equations, respectively, with nonlinear and linear stiffness and is the state vector including displacement and velocity.
In fact, there are two intersection points and of the border and the periodic orbit of primary resonance, as shown in Figure 4. Without losing generality, the neighborhood indicated by of 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 , , , and . The points in and () are accordingly governed by the solution flows of (49) and (50) in subspaces and , respectively. And and are dominated by (49) and (50), respectively. and 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 is smooth, the orbit starting from initial conditions in the Poincare section satisfies Taking the derivative of expression (53) with respect to at yields Then the derivative of (54) with respect to time gives Considering [27], the Jacobian matrix of submap satisfies the equation , where subscript denotes the time during which the orbit moves from to in subspace . Therefore, where satisfying the initial condition of is the fundamental solution matrix of (49). Considering that (49) is nonlinear with respect to state variables, it is difficult to obtain analytically. An alternative approach is a numerical method like the Runge-Kutta scheme to find the solution of the initial-value problem with Similarly, the Jacobian matrix of the other submap may be calculated in the same way.
Having established the Jacobian matrices of submaps and , the next step should be to find the Jacobian matrixes of and . Actually, is a local map representing the switch transform of state variables before and after the orbit crossing the border. And is not differentiable anywhere due to the discontinuity at the break point. In fact, (named as saltation matrix or switch matrix) symbolizes the relationship of variation to , which denote the variations between the undisturbed orbit and disturbed near the border, as can be seen in Figure 5 [27, 28]. In retrospect, and 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 reaches the border to a turn, but the disturbed orbit lags behind by a short distance to the border.
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 . Especially, the original orbit and disturbed orbit will return to the section at the same time and then set forth again, because both orbits possess the same periodicity and thus actually is an identity matrix and can be omitted. Finally,
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: , , , , , , and . 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 rad/s, the approximate analysis method in Section 3.1 finds , , , , , and . 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 , which proves that the periodic orbit is stable. This conclusion could be confirmed by the result of 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 , , , , , , , and .The procedure discussed in Section 3.2 gives , , , , , and . 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.
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.
(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, Figure 11 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.
(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 and 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.
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This work was supported by the Natural Science Foundation of China under Grant no. 11272145, Funding of Jiangsu Innovation Program for Graduate Education under Grant no. CXLX12_0135, and the Fundamental Research Funds for the Central Universities.