Numerical simulation of the soft contact dynamics of an impacting bilinear oscillator

Systems constituted by moving components that make intermittent contacts with each other can be modelled by a system of ordinary differential equations containing piecewise linear terms. We consider a soft impact bilinear oscillator for which we obtain bifurcation diagrams, Lyapunov coefﬁcients, return maps and phase portraits of the response. Besides Lyapunov coefﬁcients diagrams, bifurcation diagrams are represented in terms of both non-dimensional time instants of contact (when the mass impacts the obstacle) and of por-tions of contact duration (the percentage-time interval when the material point is inside the obstacle) vs. non-dimensional external force frequency (or amplitude). The second kind of diagrams is needed because the contact duration (or the complementary ﬂight time duration) are quantities that can easily be measured in an experiment aiming at conﬁrming the validity of the present model. Lyapunov coefﬁcients are evaluated converting the piece-wise linear system of ordinary differential equations into a map, the so-called impact map, where time and velocity corresponding to a given impact are evaluated as functions of time and velocity corresponding to the previous impact. Thus, the usual methods related to this last map are used. The trajectories are represented in terms of return maps (all points in the time-velocity plane involved in the impact events) and phase portraits (the trajec-tory-itself in the displacement-velocity plane). In the bifurcation diagrams, transition between different responses is evidenced and a perfect correlation between chaotic (peri-odic) attractors and positive (negative) values of the maximum Lyapunov coefﬁcient is found.


Introduction
Significant research efforts have been spent in the theory and application of nonlinear dynamics for non-smooth systems [1][2][3][4][5][6]. Among the wide range of nonlinear dynamical systems, Piecewise Smooth Systems (PSS) play an important role and can be classified as continuous or discontinuous PSS [5,7]. The most simple discontinuous PSS is the impact oscillator studied since 1958 for the case of an electronic bell [8] and also investigated in recent theoretical and numerical works, see e.g. [9][10][11][12][13][14]. The simplest continuous PSS is the piecewise (or bilinear) oscillator, that is the object of an important work by Shaw and Holmes [15] and is also recently analyzed in many papers, see e.g. [16,17].
Both hard and soft impacts have been considered. The exact solution of the fundamental periodic motion of a simple mechanical system with hard impact, attached to a sinusoidally excited primary mass of the system, was derived analytically by Blazejczyk-Okolewska and Peterka [18]. Blazejczyk-Okolewska et al. [19] determined the regions of periodic motions with impacts and the stability of periodic solutions of a two-degree-of-freedom mechanical system; impacts between the mass and the rigid basis were described by a coefficient of restitution. Peterka and Tondl [20] dealt with the analysis, via numerical simulations, of the motion of one degree-of-freedom mechanical system with soft impact; the restoration force is characterized by a triangle hysteresis loop and the external excitation is an harmonic force. Subharmonic impact motions are characterized in a subregion of the plane of dimensionless excitation frequency and static clearance, showing the regions of different regimes of impact motions.
Recently, various types of soft impacting systems have been addressed by Ma et al. [21] via experimental and/or numerical approaches, in the framework of a two-dimensional mapping analysis of the relevant dynamics. Among authors using sophisticated analytical mappings in the study of piecewise systems we mention Luo [22] and Pavlovskaia and Wiercigroch [23].
In this work, we also deal with a soft impact bilinear oscillator which, besides its own interest, also aims at representing a Single Degree Of Freedom (SDOF) model of the first mode of vibration of an impacted cantilever beam of uniform mass, experimentally studied in [11]. The target is to find the conditions for which non-trivial impacting chaotic and periodic motions do occur. To this aim, we use impact maps, whose limitations for studying grazing and low-velocity impacts do not come into play. In order to grossly distinguish between a periodic and a chaotic trajectory, we use Lyapunov coefficients. Lyapunov numbers (exponents or coefficients) measure the average divergence of nearby trajectories. A chaotic system is generally defined under the condition that the associated largest Lyapunov coefficient is positive [24].
For dynamical systems described by smooth differential equations and for discrete maps, the calculation of Lyapunov coefficients is well developed [25]. For non-smooth dynamical systems we have the work of Müller [26] that generalizes classical techniques for smooth dynamical systems, but also the papers of Stefanski [27], who uses chaos synchronization, and Galvanetto [28], who uses the definition of a smooth transcendental map. This map is defined in such a way the event of a certain discontinuity in the solution of the non-smooth dynamical system is given in terms of the foregoing discontinuity and it is also used by de Souza and Caldas [13]. In [13], the authors apply the usual method developed in [25] for discrete maps to this new transcendental map. They apply this idea to the case of an impact oscillator and to the impact pair system. A novelty of the present work lies in the application of the technique used by de Souza and Caldas to the case of an impacting bilinear oscillator. We remark that the tricky point of this application is that the transcendental maps to be defined are two: one map for each smooth region of the problem, namely the contact and the flight regions. Moreover, we study the stability of the trajectories by using the Jacobians in numerical form calculated at the time instants of the attachment and detachment of the material point at the obstacle, as well as the Lyapunov coefficients.
This means that it is not possible to analyze the stability of a period one point through the evaluation of its eigenvalues, see [14], without an empirical assumption on the periods of such trajectories, see e.g. [15]. This is the reason why we will evaluate numerically the stability of the trajectories by using a proper form of the Jacobian.
In order to pursue the analogy with an impacted cantilever beam, not only the springs have two different rigidities in the two smooth regions of the problem, but also the dampers have different attenuation coefficients. The spring and the damper of the system simulate the cantilever beam; the spring and the damper modelling the obstacle have larger values of rigidity (in order to simulate the hardness of the obstacle) and attenuation coefficient (in order to simulate the energy dispersion during the impact event). We remark that, in the impact oscillator model, the hardness of the obstacle and the energy dispersion of the impact event are modelled, respectively, by the fact that the latter occurs instantaneously and by a single restitution coefficient lower than one.
The paper is organized as follows. In Section 2 we present the dimensional and non-dimensional sets of piecewise ordinary differential equations that we treat in the rest of the paper. In Section 3 we describe the iterative solution method and derive the transcendental maps and their Jacobians. In Section 4 the Lyapunov coefficients are described along with the numerical method for their computation. In Section 5 we present a detailed parametric analysis of system non-trivial impact-response through return maps and bifurcation diagrams, and characterize some relevant transition scenarios. Our conclusions end the paper in Section 6.

Dimensional equations
The piecewise linear oscillator considered in this paper is a mass-spring-damper system governed by the following ordinary differential equation, where x is the displacement, from the unstressed configuration of the spring with rigidity k s , of a material point with mass m at time t, the dot denotes the derivative with respect to time, F 0 sinðx tÞ is the external sinusoidal force; cðxÞ is the attenuation coefficient and kðxÞ is the opposite of the force exerted by the system springs. In order to simulate the occurrence of the impact event, cðxÞ and kðxÞ are two piecewise functions of x defined as follows, where d is the amplitude of the gap, that is the distance of the point, from the obstacle, in the unstressed configuration of the spring with rigidity k s . The pedix s refers to the system composed of the spring with rigidity k s and the damper with attenuation coefficient c s . The pedix o refers to the obstacle; the obstacle is simulated by the spring with rigidity k o , whose unstressed configuration is shifted by the gap d from the unstressed configuration of the spring with rigidity k s , and by the damper with attenuation coefficient c o . Both the spring with rigidity k o and the damper with attenuation coefficient c o exert forces, respectively Àk o ðx À dÞ and Àc o _ x, only if x > d. Generally, see e.g. [15], the attenuation coefficient related to the obstacle is assumed to be zero, c o ¼ 0. However, during the impact, the energy dissipation is larger than with no impact. Thus, the assumption will be considered in this paper. A graphical scheme of the model described by the Eqs. (1)-(3) is given in Fig. 1. For the piecewise linear ordinary differential Eqs. (1)-(3) an analytical solution is available in the regions where xd (where (1)-(3) are smooth) and once initial conditions are given. The tricky point is the shifting from x < d to x > d and viceversa, whose time instants are herein evaluated numerically, consistent with the overall framework of the paper, although analytical solutions could also be pursued under proper conditions, see e.g. [18,29].

Non-dimensional equations
Characteristic dimensions of this piecewise linear system are the gap d, the mass m of the point and the rigidity k s of the system's spring. With these three values we define all the non-dimensional quantities we use in the paper, that will be denoted with a tilde over the analogous dimensional symbol, i.e., so that the time derivative operators and the functions defined in (2) and (3)  x þk o ðx À 1Þ;x P 1; ( ð8Þ cðxÞ ¼c s ;x < 1; c s þc o ;x P 1:  A similar non-dimensional-procedure is also used in [12]. From now on, for the sake of simplicity and except when differently mentioned, we will suppress the tilde over the symbols, and the dimensionless time derivative operators will be denoted like the usual time derivative operators, i.e., by dots over the functions to be derived; nevertheless, all the quantities will be considered non-dimensional according to the rules of this subsection.

Poincaré sections, return map and phase portrait
For the 2nd order non-autonomous ordinary differential Eq. (1) the phase space P is three-dimensional ðt; x; _ xÞ 2 P. A natural 2D section of this space is made at x ¼ 1 where the discontinuity occurs, see Eqs. (7)- (9). The Poincaré section R r considered in this paper is defined as follows, so that a point of the Poincaré section R r is represented by the couple ðt m ; _ x m Þ 2 R r such that, Let us remark that the modular operation is done without lack of generality because the time t appears explicitly in the governing equation only in the function sin xt, on the right-hand side of (7), that is periodic of the same period 2p=x considered in the modular operation. This modular operation allows us to represent the points of the Poincaré section R r in the compact graph R r which is called the return map.
Another place to slice the phase space P is at t ¼ t c > t 0 (with t c constant). In this case, the Poincaré section R tc is defined as follows, For every dimensionless time constant t c , we have one point in the plane ðx; _ xÞ. If one uses t c as a curve parameter the graph is generally called the phase portrait.

Preliminaries to the scheme of the iterative solution
If the initial conditions are chosen in such a way the point does not touch the obstacle ðx -1Þ, than the solution of the problem is trivial before the first impact. Otherwise, with no lack of generality we choose the initial condition on the gap, i.e., at the initial dimensionless time t ¼ t 0 , we have x ¼ 1. If the initial dimensionless velocity is positive _ x 0 > 0, then the ordinary differential equation to be solved is obtained from (7)-(9) as If the initial velocity is negative _ x 0 < 0, then the ordinary differential equation to be solved is obtained from (7)- (9) as Differential Eqs. (13) and (14) or (15) and (16) can be solved analytically.
The point 0 of the return map R r is given by the initial condition ðt 0 ; _ x 0 Þ. The point 1, ðt 1 ; _ x 1 Þ, is deduced by the solution xðtÞ of (13) and (14) (if _ x 0 > 0) or of (15) and (16) (if _ x 0 < 0) and it is defined by the following proposition, The point 2 ðt 2 ; _ x 2 Þ of the return map R r is derived in the same way. Let us remark that the sequence ð _ . . .Þ of dimensionless velocities on the return map R r has alternation of signs because if for the ith point of R r the material point penetrates the obstacle, then for the ði þ 1Þth point of R r the material point comes out of the obstacle and vice versa. The possibility of having null initial velocity is excluded by (17); in other words the phenomenon of grazing is included in the phase of flight motion, because we are interested in treating non-trivial cases of obstacle penetration. Moreover, since the sequence ðt 0 ; t 1 ; t 2 ; . . .Þ of dimensionless times on the Poincaré section R r is strictly increasing, in the return map R r defined in (11) this property is bypassed because of the modular operation. The formal presentation of the numerical simulations made in the paper is given in the next subsection.

Scheme of the iterative solution
Without lack of generality, let us assume that the point ð2iÞ ðt 2i ; _ x 2i Þ on R r is such that _ x 2i > 0. In order to evaluate the point ð2i þ 1Þ ðt 2iþ1 ; _ x 2iþ1 Þ on R r , the differential Eq. (13) with initial conditions (14) gives a solution expressed with the function x o , (read o as obstacle), defined as follows, , the dimensionless velocity of the point and t 2iþ1 is the first dimensionless time instant (after t ¼ t 2i ) when the dimensionless displacement xðtÞ reaches again the value of the gap The analytical expressions of the functions x o and t o are complicated but they can be determined via a clear pattern (see e.g. Section 3.2). Therefore, the value of the dimensionless velocity t o at dimensionless time so that the point ð2i þ 1Þðt 2iþ1 ; _ x 2iþ1 Þ on R r is evaluated. In order to evaluate the next point ð2i þ 2Þðt 2iþ2 ; _ x 2iþ2 Þ on R r , the differential Eq. (15) with initial conditions (16) gives a solution expressed with the function x s , (read s as system), defined as follows, where t s ðt; t 2iþ1 ; _ x 2iþ1 Þ is, inside the dimensionless time interval ½t 2iþ1 ; t 2iþ2 , the dimensionless velocity of the point and t 2iþ2 is the first dimensionless time instant (after t ¼ t 2iþ1 ) when the dimensionless displacement xðtÞ reaches again the value of the gap x ¼ 1, The analytical expression of the functions x s and t s are complicated but determinable through a clear pattern (see e.g. the Subsection 3.2). The value of the dimensionless velocity t s at dimensionless time t ¼ t 2iþ2 is denoted by _ so that the point ð2i þ 2Þðt 2iþ2 ; _ x 2iþ2 Þ on R r is evaluated.
The process described above is iterated for every i. We start with i ¼ 0 and from the point 0 ðt 0 ; _ x 0 Þ on R r we evaluate the points 1, 2 and so on.

The return map and its Jacobian
The procedure explained in the previous subsection is used to define the following return map, where P is an operator that takes a point ðt 2i ; _ x 2i Þ on R r (with positive dimensionless velocity) and gives back the next point ðt 2iþ2 ; _ x 2iþ2 Þ on R r (with positive dimensionless velocity). In order to evaluate the form of the return map P, the first step is to obtain the couple ðt 2iþ1 ; _ x 2iþ1 Þ as a function of ðt 2i ; _ x 2i Þ, where the functions F o and G o are deduced from the functions x o and t o in Eq. (18). The second step is the evaluation of the couple ðt 2iþ2 ; _ x 2iþ2 Þ as a function of ðt 2iþ1 ; _ x 2iþ1 Þ, where the functions F s and G s are deduced from the functions x s and t s in Eq. (23). Thus, the return map is completely defined as follows, where the functions F and G are the two components of the return map P and are self-defined in (33) and (34) in terms of the functions F o ; F s ; G o and G s . It is interesting to evaluate the Jacobians J o ; J s and J of the maps defined in (29,30), (31,32) and (33, 34), respectively, An analytical and explicit expression of the maps F o and F s can not be found because the conditions from (20) and (25), are generally not invertible in our case. In fact, our investigation is not limited to find periodic solutions of motion, see Blazejczyk-Okolewska et al. [19]. On the other hand, the derivatives of (36) can be analytically computed, once the form of the functions F o ðt 2i ; _ x 2i Þ and F s ðt 2iþ1 ; _ x 2iþ1 Þ are used in (36), Differentiating the expressions in (37) with respect to t 2i and _ x 2i and with respect to t 2iþ1 and _ x 2iþ1 , respectively, we obtain, that can be inverted to obtain the first rows of the Jacobians J o and J s , On the other hand, an analytical and explicit expression for the maps G o and G s are obtained evaluating the velocity functions t o and t s , respectively at t ¼ t 2iþ1 and at t ¼ t 2iþ2 , i.e., so that the functions G o and G s are defined through the functions t o ; t s ; F o and F s , and the second rows of the Jacobians J o and J s are, @G s @t 2iþ1 ¼ @t s @t 2iþ2 @F s @t 2iþ1 þ @t s @t 2iþ1 ¼ À @ts @t 2iþ2 @xs @t 2iþ1 @xs @t 2iþ2 þ @t s @t 2iþ1 ; ð52Þ The expressions of the Jacobians J o and J s are tedious and long to be written explicitly but they are used in the numerical code. For the chain derivative rule, the expressions of the components of the Jacobian J are, so that the Jacobian matrices defined in (35) can be arranged as follows, where the symbol Á denotes the classic rows by columns product. In order to achieve the form of the Jacobian matrix as a function of t 2i and of _ x 2i , the evaluation of the matrix J s should be done at the same point. To do this, Eqs. (29) and (30) must be used, However, as we have already pointed out, an analytical expression of the function F o is not available in its general form. Thus, the Jacobian matrix J 2i ¼ Jðt 2i ; _ x 2i Þ can not be calculated analytically. This means that it is not possible to analyze the stability of a period-one point through the evaluation of its eigenvalues, see [14], without an empirical assumption on the periods of such trajectories, see e.g. [15], in a soft contact system. This is the reason why, on the other hand, we evaluate numerically the stability of the trajectories by using the form of the Jacobian J 2i given in (54).

Error propagation in the return map
If one assumes to have an error e 2i ¼ ðdt 2i ; d _ x 2i Þ > in the evaluation of the point 2i ðt 2i ; _ x 2i Þ on R r , then the error e 2iþ2 ¼ ðdt 2iþ2 ; d _ x 2iþ2 Þ > in the evaluation of the point 2i þ 2ðt 2iþ2 ; _ x 2iþ2 Þ > on R r can be computed through a Taylor's series expansion, Thus, for an uncertainty e 0 ¼ ðdt 0 ; d _ x 0 Þ > in the evaluation of the initial point 0 ðt 0 ; _

Lyapunov exponents
We describe in this section the computation algorithm to calculate the Lyapunov exponents, that are defined as follows, where K N h is the hth eigenvalue of the matrix M 2N defined in (58).
Many authors [13,24,25] show that for chaotic attractors the components of M 2N , even for few iterations, become so large that numerical errors are huge and the evaluation of Lyapunov exponents must be approximated. The procedure is illustrated in [25] and schematically presented here; the Jacobian matrix J 0 is decomposed into the product of an orthogonal matrix Q 0 and an upper triangular matrix T 0 with non-negative diagonal elements. If J 0 is invertible, the decomposition is unique and we have, The other Jacobian matrices J 2i with i ¼ 1; 2; . . . ; N À 2 are decomposed, iteratively, as follows, in such a way we have, Following [25], the assumption Q 2NÀ2 equal to the identity matrix is done and (59) can be replaced by the following computationally friendly approximation, where a 2j h is the hth non-negative diagonal element of the upper triangular matrix T 2j . Let us note that with the approximation Q 2NÀ2 ¼ I, Eq. (63) is derived considering the fact that the product of triangular matrices T j is a triangular matrix (the diagonal elements of which are the products of the diagonal elements of the triangular matrices T j ), that the eigenvalues of a triangular matrix correspond to the elements of the diagonal and that the logarithm of a product is the sum of logarithms.

Parametric analysis
To perform a parametric analysis, we have to set a certain number of parameters. To do this we refer to an impacted cantilever beam with uniform mass distribution in the background. The identification with the present Single Degree Of Freedom (SDOF) model can be done assuming that the tip displacement of the beam under a static force is the same of the displacement x corresponding to the same applied static force, and that the resonance frequency of the SDOF model and the first eigenfrequency x I of the beam model are equal, k s ¼ 3EI where E is the Young modulus of the beam material, I is the moment of inertia of its section with respect to the principal axis orthogonal to the plane of motion, L is its length and the dimensional notation is resumed. A specific case is shown in Table 1 and the non-dimensional value of the system's attenuation coefficient c s is from (5) 4 , Attenuation coefficient c s is chosen in such a way the oscillation amplitude is attenuated at the rate of 10% every 10 cycles. The remaining dimensionless parameters are the rigidity k o and the attenuation coefficient c o of the obstacle, the amplitude F 0 and the frequency x of the external force. In this bilinear oscillator (a soft impact model) the restitution coefficient r is not assigned a priori as in a hard-impact model. Indeed, in this soft impact model the modulus of the ratio of the velocity t a after the impact and t b before the impact can be numerically computed. In impact oscillator (hard-impact) models the interval of time Dt in during which the mass is inside the obstacle is fixed to be zero and the interval of time Dt out between successive impacts is called flight time. In this (bilinear oscillator) model with soft impacting obstacle t b ; t a ; r ¼ Àt a =t b ; Dt in and Dt out are quantities that can be computed in the numerical simulation. In the experiments shown in [11] on an impacted cantilever beam, it is shown that the restitution coefficient r and the intervals Dt in and Dt out assume approximately the following values, Dt in Dt out ' ð1-2%Þ: ð66Þ Table 1 Young modulus, moment of inertia, length and first eigenfrequency of the beam used in the experiments in [11] (dimensional values). From (64) m and ks are evaluated.
that will be used in the simulations. In Fig. 2, return maps and phase portraits are shown in the first and in the second row, respectively. In the six graphics of Fig. 2, numerical values are the following: k o ¼ 10 3 ; k s ¼ 1; c o ¼ 6:2; c s ¼ 6:2 Â 10 À4 ; F 0 ¼ 1:4881 and m ¼ d ¼ 1. Three different sample cases of response are reported for the frequency values x ¼ 2:211 (1st column) , x ¼ 2:176 (2nd column) , and x ¼ 3:750 (3rd column). In the first column ðx ¼ 2:211Þ an attractor of period eight is seen, so that the 16 points represented in the return map (8 for dimensionless time instants of contact and dimensionless velocities at the beginning of the contact and 8 at the end of the contact) are the points of the attractor and the closed orbits in the phase portrait are exactly eight. In the second column ðx ¼ 2:176Þ the attractor is strange, it is composed in the return map by an infinite number of points and in the phase portrait by an infinite number of closed orbits. In the third column ðx ¼ 3:750Þ a perfect stationary (one period) case is presented: the attractor is made of 2 points represented in the return map (1 at the beginning and 1 at the end of the contact) and the unique closed orbit in the phase portrait is shown.
The bifurcation analysis is done by considering the non-dimensional external force characteristics x and F 0 as control parameters. In Figs. 3 and 4 the non-dimensional external force amplitude is fixed to be F 0 ¼ 1:4881, whereas in Figs. 5 and 6 the non-dimensional external force frequency is fixed to be x ¼ 1, both corresponding, e.g., to the dimensional values of Table 1. The details of the results in Figs. 3-6 are expressed in the following comments.
In Fig. 3 the parameterization is done with respect to the non-dimensional external force frequency x. For each value of x, in the first row, the two Lyapunov exponents are given. In the second row, the non-dimensional time instants of contact when the mass recovers the gap and impacts the obstacle at x ¼ 1 are reported. In the third row the portion of contact duration (the time interval Dt in when the point is inside the obstacle) is calculated in percentage, see Eq. (66) 2 , of the total time of contact duration ðDt in Þ plus the successive flight time duration (the time interval Dt out when the point is outside the obstacle). In the second row, two different symbols are used to indicate the inward () and the outward (h) directions with respect to the obstacle, see also the enlargements in Fig. 4. The time instants of contact are modulated according to the period 2p=x of the external force. For the sake of clearness, the first and second column display two different ranges of x; in the first column, we show the range between x ¼ 0 and x ¼ 2:2; in the second column, the range between x ¼ 2 and x ¼ 5. In Fig. 4 we enlarge the boxed parts of Fig. 3 so that the non-dimensional frequency range is between x ¼ 3:32 and x ¼ 3:38. We provide detailed analysis of what happens in this very narrow range. Such analysis has limited importance for engineers designing impacting systems; however, it is of great interest for understanding transition scenarios between various kinds of motion. decreasing values of x) we have continuous and discontinuous bifurcations, see e.g. [5,30], as well as chaotic behavior, evidenced by the presence of a positive maximum Lyapunov exponent, and very short flight phases. Also for x 2 ð3:76; 4:384Þ the point remains confined in the regular region x < 1 where no contact occurs (flight phase). For x 2 ½3:323; 3:76 the point has contact with the obstacle; in the right part of this range ðx 2 ½3:6; 3:76Þ we have a period 1 orbit ending at x ¼ 3:6 with a continuous bifurcation. One of the eigenvalues of the Jacobian is equal to À1, so that the bifurcation is a period doubling, see e.g. [31]. Thus, for x 2 ½3:39; 3:6Þ we have a period 2 orbit. For x 2 ½3:372; 3:39Þ, see Fig. 4, there is a sequence of discontinuous bifurcations giving rise to alternating period 2 and period 3 orbits; given their dense pattern in a very narrow frequency, it is considered not worth to characterize these bifurcations through the Jacobian of the map, along the lines of [22].    Fig. 4, we have a continuous bifurcation from period 2 orbits to period 4 orbits. In the range x 2 ½3:354; 3:371Þ we have a period 4 orbit and in the range x 2 ½3:343; 3:354Þ again a sequence of discontinuous bifurcations characterized by an alternation of higher periodicity orbits in a very tiny region. However, for x ¼ 3:344 the attractor is strange and the relevant maximum Lyapunov exponent is positive. Even in the range x 2 ½3:323; 3:342 the attractor is strange and the maximum Lyapunov exponent is positive. However, for x ¼ 3:326, the attractor is not strange and the maximum Lyapunov exponent is negative. For x 2 ð2:77; 3:323Þ, see Fig. 3 (second column), the point remains confined in one regular region, i.e., x < 1, where no contact occurs (flight phase). In the range x 2 ½2:152; 2:77 the point has contact with the obstacle and the qualitative behavior is analogous to that in the range x 2 ½3:323; 3:76. In the range x 2 ½1:77; 2:152, see x ¼ 2:11. In the range x 2 ½0:81; 1:77 there is a period 1 orbit with a discontinuous bifurcation occurring at x ¼ 0:81. Then, for x 2 ½0:51; 0:81 we have a period 2 orbit and, finally, for x < 0:51, an alternation of continuous and discontinuous bifurcations as well as of periodic and strange attractors. Overall, looking at the whole decreasing frequency range, a structured sequence of non-impacting and impacting regions is observed. The former ones are followed by regions of periodic impacting motions, which bifurcate to chaos through a sequence of continuous (period doubling) and discontinuous bifurcations, with the chaotic region ending up nearly abruptly to non-impacting dynamics.
In Figs. 5 and 6 the parameterization is done with respect to the non-dimensional external force amplitude F 0 . Also in this case the agreement between regions of positive (negative) largest Lyapunov exponent and chaotic (periodic) impact motions can be observed. Two large ranges of non-dimensional driving force amplitude, one where the response is stationary (between F 0 ¼ 0:2 and F 0 ¼ 6:13, see Fig. 5) and one where is chaotic (between F 0 ¼ 6:3 and F 0 ¼ 40:88, see Fig. 5) are exhibited, as well as two very small in-between ranges with the same characteristics (between F 0 ¼ 6:16 and F 0 ¼ 6:29 for the stationary case and between F 0 ¼ 6:13 and F 0 ¼ 6:16 for the chaotic case, see the 1st column of Fig. 6). The abrupt passage from periodic to robust chaotic orbit is evidenced in the first column pictures of Fig. 6, where the range between F 0 ¼ 6 and F 0 ¼ 6:4 is shown. In the second column pictures of Fig. 6, another chaotic range between F 0 ¼ 40 and F 0 ¼ 42 is shown; herein, the solution is confined in the regular region x < 1 (detachment) for the non-dimensional driving force amplitude between F 0 ¼ 40:88 and F 0 ¼ 41:08.
In view of the identification with the beam impacting problem in the background, we remark that the absolute value of the gap can also be considered as a measure of the damage of a certain part of the cantilever beam. Thus, it is clear the importance to perform a parametric analysis also with respect to d. Actually, in the non-dimensional context of this section, such a parameterization is indirectly performed via the F 0 -parameterization because of Eq. (5) 6 .

Conclusions
For a sinusoidally forced bilinear oscillator (with the discontinuity displaced by a gap from the unstressed configuration of one spring) we implement a numerical iterative method to derive the solution and discuss some relevant aspects. The dynamics of the soft impact oscillator, which is a discontinuous PSS, is investigated based on a return map obtained by slicing the three-dimensional phase space at the discontinuity. The definition of this transcendental map is not explicit, but an explicit analytical form of its Jacobian is given along with its evaluation on the computed trajectory. We use the Jacobians to evaluate the Lyapunov coefficients by the method exposed in [13]. We show phase portraits and return maps for three different response samples: when the orbits are infinite in number (chaotic behavior and a positive largest Lyapunov coefficient), when the orbits are finite in number (high-periodic response with a negative Lyapunov coefficient) and when the orbit is unique (period-one response corresponding again to a negative Lyapunov coefficient).
Besides classical bifurcation diagrams in terms of Lyapunov coefficients, more characterizing diagrams are given in terms of the time instants of contact when the mass recovers the gap and impacts the obstacle, as well as of the percentage-portions of contact duration. The former are shown to grossly distinguish between regions of periodic and chaotic response. The latter allow us to evaluate bifurcation points, discuss routes to chaos and detect periodic orbits, while better understanding the mechanical features of the observed (in/out of the soft contact) dynamics. Overall, a very complex scenario has been highlighted, with a rich structured pattern that includes continuous and discontinuous bifurcations, flip bifurcations and alternating regular and chaotic regimes. The importance of the second kind of diagrams is also somehow related to the possibility of an experimental investigation on an impacted cantilever beam, where the measure of contact duration and of flight time is easy to obtain, as shown in [11]. Indeed, the considered bilinear oscillator can also represent a simple SDOF model for effectively describing the impacting behavior of a cantilevered beam with uniform mass distribution. It can be useful (i) to explore a wide range of problem parameters with low computational efforts, (ii) to capture interesting features of the dynamic response, and (iii) to get hints for further investigations of the beam via richer reduced order models [19] or 1D FEM models accounting for intermediate loading positions along the beam axis, for coupled effects of other contacting supports, and for the contribution of higher modes. More generally, interpreting the clearance extent between mass and obstacle as structural damage and analyzing the nonlinear characteristics of the forced response of the system could provide valuable means for damage detection via vibration-based methods.