J P Cusumano

B.-Y Bai

Period-infinity Periodic Motions, Chaos, and Spatial Coherence in a 10 Degree of Freedom Impact Oscillator

The numerical study of a 10 degree of freedom impact oscillator is presented. The mathematical model consists of a linear array of masses in which each mass is connected to its nearest neighbor by identical springs and dashpots. One end mass is attached to a rigid support while the other is free to impact a sinusoidally vibrating rigid table. Bifurcation diagrams based on the impact Poincare map are obtained over the entire range of natural frequencies, taking the table frequency as the bifurcation parameter.

The diagrams reveal many chaotic bands as well as a wide variety of period-(n,m) (P,l,,,) orbits, where n is the period with respect to the impact PoincarC map and m is the number of table periods.

Perhaps the most interesting solutions are P,,,,, solutions (where m is finite). These arise from sticking events in which the time between impacts approaches zero in finite time. The convergence properties of the impact time sequence during a sticking event is shown to be controlled by the coefficient of restitution of the impact law. The spatial complexity of the motions is studied by finding proper orthogonal modes for various steady states, and the number of excited degrees of freedom is estimated by finding the number of modes needed to exceed 99% of the signal power. These results are compared to the predictions of fractal dimension theory.

INTRODUCTION

Impact oscillators

have furnished a useful mechanical paradigm for discrete-time nonlinear dynamical systems, and have been the subject of many recent investigations, both theoretical [l-6], and experimental [7-l 11. In engineering systems, impact is sometimes an essential part of the design, as in impact printers [12]; but, more often, it is the source of unwanted wear and noise, as in the dynamics of structures with joint clearances, or in gear systems. Whereas most previous studies of the nonlinear dynamics of impact oscillators have involved systems with, either literally or effectively, one degree of freedom, flexibility effects should be expected to play an important role in applications. Furthermore, from the elementary theory of impact, one anticipates that many degrees of freedom will be excited by impacts on flexible structures.

Thus, flexible impact oscillators appear to be an ideal test bed with which to explore ideas of dimensionality and spatial complexity in the dynamics of solids and structures.

In this paper, we present the results of a numerical study of a 10 degree of freedom impact oscillator.

The system is a linear array of masses in which each mass is connected to its nearest neighbor by identical linear springs and dashpots and in which one end mass is free to impact a sinusoidally vibrating rigid table. Since all of the nonlinearity is contained in the impact boundary condition, analytical solutions are available for motions between impacts. These solutions are used to construct the implicit 20-dimensiorial Poincark map defined by the impact, and a computer code based on this map is used to obtain our results.

In Section 2, we describe the physical system, its continuous-time mathematical model, and the associated impact-based Poincark map. In Section 3, details of the computer 515 1 algorithm used in our numerical experiments are presented, and the properties of sticking events are described: it is shown that the impact time sequence during a sticking event is controlled by the coefficient of restitution of the impact law and that sticking occurs in finite time. A modified algorithm is described which allows the simulations to proceed across sticking events. The bifurcation structure of the system is presented in Section 4, and nature of the various solutions are described with the help of Poincare sections, including chaotic as well as a wide variety of period-(n.m) (P,,,,,)) solutions, where n is the period with respect to the impact Poincare map and m is the number of table periods. In particular, it is shown that sticking events can result in P,,,,, solutions (where vz is finite). In Section 5, we give a brief introduction to the proper orthogonal decomposition method for the study of spatial complexity, and we describe its application to our problem. The results of the proper orthogonal decomposition analysis are then presented, along with fractal dimension estimates for various steady states, and results are compared to the predictions of dimension theory. Section 6 is reserved for a summary and discussion of our results.

THE MATHEMATICAL MODEL AND METHOD OF SOLUTION

The object of this study is the system shown in Fig. 1. It consists of 10 identical masses of mass m connected to nearest neighbors by identical linear springs and dashpots with constants k, and c, respectively. One end mass is connected to a rigid wall, again with a linear spring and dashpot.

whereas the other end mass is free to impact on a sinusoidally moving rigid table. The entire system is constrained to permit only horizontal motion, with coordinates X, through x,,, defined in an obvious manner as in the figure. Between impacts, the system is completely linear, and thus can be represented by a dimensionless system of equations in the form S + ;fAi + Ax = 0 [START_REF] Holmes | The dynamics of repeated impacts with a sinusoidally vibrating table[END_REF] where xT = (x,, x2,. . ., xl,,), y = cq(km) is a dimensionless damping parameter, and the overdot denotes the derivative with respect to the dimensionless time t, which is related to the dimensional time i by t = (v'( k/m));.

A is a 10 x 10 banded matrix given by: A=

2 -1 0 -1 2 -1 0 -1 2 0 0 0 0 0 0 0 0 (I 0 0 o- 0 0 0 0 0 0 2 -i 0 -1 2 -1 0 - 1 1
(2)

Using the modal matrix @ (normalized so that (DT<D = I). we transform the system to modal coordinates qT = (q,, q2, . , c/,,,) by ., lo), [START_REF] Whiston | Global cfvnamics of a vibro-impacting linear oscillator[END_REF] in which the modal damping ratios <, are, as usual, defined by C, = (yo;)/2. The solution to the system of equations ( 5) is then simply 9, = ee:rCU! (a, cos b,t + b, sin CJ$), [START_REF] Whiston | Singularltles in vibro-impact dynamics[END_REF] where I& = oll/(l -cf) is the damped natural frequency, and the constants a, and bi are determined from the initial data at time to as a, = ql(fo); and b, = 4Xto) + C,w,q,(to) a, Given the modal solutions ( 6), the solution in spatial coordinates is then obtained from the transformation (3). (The natural frequencies for the system between impacts are presented in Appendix A. We defer display of the normal modes until the discussion of spatial coherence in Section 5.) The system evolves according to the solutions (6) until an impact occurs. That is, until

x,o(tJ = XT(&). ( 8
)
where t,, is the time of the nth impact and the table motion is given by xr( t) = A sin ot + A,,. Equation ( 8) defines a scalar-valued function of time, the first root of which gives the impact time. Employing a simple impact theory in which the impact is assumed to occur instantaneously, just after the nth impact one has

x,(tl:) = x,(t,) = x,(t,,) (i = 1, 2, . . .) 10) (9. where u, denotes the velocity of the ith mass. and plus and minus superscripts indicate the right and left hand limits as t + t,,, respectively.

In the impact law (9.3), cy is the coefficient of restitution, and we use the fact that +(tz) = QtJ. Relations (9) give new initial conditions at time t,,, which are then used to generate the solution for the system up to the next impact, via equations (3), [START_REF] Whiston | Singularltles in vibro-impact dynamics[END_REF], and (7). Thus, the dynamics of the system is completely determined by the state just after each impact, and one need only consider the 20-dimensional Poincare map defined by the impact condition, Y: IX + X, where the cross section C is given by

X = {(x, v, t) E R'" x R'" x S'/C,(,~X = XT, e^,,*v > UT}, ( 10
)
and where it is understood that I is taken modulo 2n/(0 (see Fig. 2). Let {Cl}:!!, denote the standard basis, so that x. e^, is the projection onto the ith component of x. Then, using the impact conditions ( 8) and (lo), 8 can be written as:

x,,,,+i = e^, * 'Dq(t,, + At,,; x,,, v,,) (i = 1, 2, . . ., 9) (11. where At,, = A t,,(~,~, v,, , t,,) is determined by using equation ( 8) to find the time of the (n + l)th impact. Using the same notation as equations (ll), equation ( 8) becomes D(At,,) -C,,, . @q(t,, + At,,; x,!, v,,) -A sin(t,, + At,,) -A,, = 0. (12)

The functional dependance of q and cj on the state at the nth impact, (x,~, v,,), comes from the evaluation of the constants in equations ( 7) with t,, = t,,, and from the transformation (3). 12) is very efficient as long as the time between impacts is not too small. However, in certain orbits, sticking events can occur in which the At, rapidly decrease by many orders of magnitude:

the calculations drop into the round-off-induced numerical noise floor, resulting in spurious, noisy bounces which decay no further (with a At of the order of 10-h-10~7).

The problem is not simply one of limited precision, however: even with infinite precision, these motions kill the algorithm, since for all practical purposes the code can never get through the event. in a kind of numerical version of Zeno's paradox. The most natural way to compensate for this is to modify the basic algorithm to assume that values of At,, below some threshold value indicate that the tip mass has 'stuck' to the moving table (i.e. (I,,, -t1.r = 0). However, the main question is whether or not these sticking events actually occur: that is, whether or not z{At,,} < x near a sticking event. To explore this further, define the impact time ratio i.e. within the noise floor.) Hence, by the ratio test, the sum of time intervals converges, and therefore the tip mass sticks to the table in finite time. Thus, sticking events are an essential aspect of the dynamics of impact oscillators, not merely ad hoc assumptions needed to make simulations feasible.

Following these observations, the minimum time step At,,, was set at lo-": should the root-finder compute a value of At,, < Armi", the tip mass is assumed to be constrained to lie on the moving table. Thus, the impacting system of Section 2 is replaced by a simple, non-impacting, 9 degree of freedom linear oscillator driven by a moving support. The time-step error E that results from this approximation is For this work, LY = 0.85, so the absolute precision of our time-step calculations was on the order of 7 x 10p6. As we discuss in Section 5, analysis of the Poincare section data showed that the long-time cumulative accuracy of our time-step calculations was of the order of 3 x lo-". By comparison, the maximum time steps were of the order of several table periods (i.e. with an order of magnitude of around lo), thus indicating a dynamic range for our time-step calculations of 70-90 dB, comparable to experimental observations possible under ideal circumstances with 12-14 bit data acquisition. The reduced system (with the end mass stuck to the table) is easily derived by deleting the 10th row of equation (l), and setting xl0 = xr = A sin it + A": the reaction force FR between xIo and the table is found to be

FR = XT + y(iT -i-9) + (XT -x9> (15)

where x9 is determined from the solution of the constrained system. The root-finder is then used to solve for the time when FR passes through zero (decreasing), thus indicating that tip and the table lose contact and the system is again governed by equations ( 11) and (12). 5

It is perhaps not too surprising that the time-step ratio A,, converges to the coefficient of restitution.

The limiting behavior during sticking events is thus identical to that seen in the elementary problem of a ball bouncing on a stationary floor [14]: when the time between impacts is very small in our system, the table and all but the end mass are essentially stationary.

That is, the system behaves like a single degree of freedom impact oscillator during a sticking event. Mathematically, the problem of sticking arises from a degeneracy in the Poincare map 9' [equations (11) and ( 12)]. The section Z is indeed global, however, it is not closed under the action of 9: there exist initial conditions X0 = (x0, vg, to) E C such that Y(X,)

E C for all n < a, but lim'Y"(X,) 4 2 since e^ro -v,, + I+.

BIFURCATION STRUCTURE

Bifurcation diagrams were constructed for the system using the table frequency (r) as the bifurcation parameter.

Unless otherwise noted, system parameters used in the simulations were as follows: the coefficient of restitution CY = 0.85; the table amplitude A = 0.10; the table offset A0 = 0; and the damping parameter y = 0.1. The algorithm used was capable only of following stable solutions: after allowing transients to die off by skipping the first 700 impacts, data from the next 100 impacts was stored for later display (in several test cases, it was found that skipping 1400 impacts before sampling left the results unchanged, so the lower number was used). The maximum value of Aw was 0.002: again, it was found in several test cases that doubling the resolution did not significantly change the results (more discussion about the stability of our results with respect to changes in Ao can be found below). Data from the steady-state taken at o = o,,~ was used to start the computations for w = w,,+i. The entire numerical experiment was started by releasing the system from rest in the shape of the first mode with unit amplitude.

Tests done by releasing the system with modal amplitudes of 0.1 and 0.5 lead to the same results: however, given the large size of the phase space and time-consuming nature of the computations, it was not possible to search effectively for disconnected branches of solutions.

To cover the entire range of natural frequencies, we take 0.14 s o z 2.0. The resulting bifurcation diagram is shown in Fig. 4, which plots the steady-state relative tip velocity U,,I -010 -uT just after impact against W, along with labels indicating the type of response (the way in which the response type was determined is discussed below). Various types of period-(n.m) (P,,,,,) orbits are labelled, where y1 is the period with respect to the number of table impacts and m is the period with respect to the number of table periods. In addition, many chaotic bands are visible (labelled by C,,). That these labels by no means exhaust the range of responses is demonstrated in the blown-up portions of Fig. 4 which are presented in Fig. 5. Figures 4 and5 demonstrate the sensitivity of the system as w ranges over certain values: this is not just a computational artifact, since the overall structure of the bifurcation diagrams in both figures are the same, even though AW is 2.0 x 10-j in Fig. 4 and2 x lo-" for Fig. 5.

Figure 6(a) displays the same data as Fig. 4, however, following Geist and Lauderborn

[15], we use the energy E (dimensionless in our case) to characterize the system just after impact, where which can be rearranged to give the bound:

(17) AE = ,& 1
(1 + a) !%!! + 0 u,l ' ii

HI c &l + @)M, (18)
"T "T

where, for small amplitude oscillations (in which, just before impact, lo$ << ic~~l), A4 = suP(t",,I/uT).

Th us energy input is controlled by the table velocity, and a better measure of input amplitude is given by the velocity amplitude show that the orbit is period-4 with respect to the table motion. We therefore conclude that this is a P,,, orbit: this seemingly paradoxical result is due to the intrinsic degeneracy of impacting systems, and can only be modelled with the system of equations ( 11) and (12) if they are modified as in our algorithm to pass through sticking events. Other P,,,, orbits can be found in a similar manner. directly for sensitive dependance on initial conditions: results obtained with initial data with components differing by 10e8 were compared, and, as expected, the periodic orbits showed no change, but chaotic orbits diverged after about 100 impacts.

DIMENSIONALITY AND SPATIAL COHERENCE

That the dynamics of spatially continuous systems may be described by the evolution along a low-dimensional attractor is a primitive concept of the 'dynamical systems' methodology.

The development of fractal dimension theory and, in particular, of the correlation dimension algorithm [16], has allowed dimensionality to be studied in physical and numerical experiments on a variety of fluid-and solid-mechanical systems (numerous examples can be found in the book by Moon [17], and in the references therein).

From our point of view, the utility of fractal dimension estimates lies in the following: given a motion on an attractor with fractal dimension d, the number of phase space dimensions m needed to contain the attractor is bounded by attractor of dimension d will leave the attractor lying in a subspace of dimension no greater than 2[d] + 1. Specifically, Marie's theorem states that the coordinatizations for which this bound is true are residual in the space of all coordinatizations: however, while residual sets are typical (i.e. generic) in a topological sense, the results we present below demonstrate that this theorem may not be as definitive as it first appears. The application of dimension theory also suffers from several other limitations:

[d] G m S 2[d] + 1, (19)
there are no a priori error estimates for given by the dyadic product: R(s, s') = (x(s, t)x(s', t)).

where x is an appropriate field variable in space s and time t (such as displacement in solid mechanics, or velocity in fluid mechanics), and the brackets ( ) denote the time average. In our case. the spatial correlation is the tensor given by: R = (x(t)x(t)) = ( ,( ) ,( ) X f X t e^,e^,) = (X,(t)X,(f))e^,e^, = R,e^,e^,, where q(t) is the displacement of the ith mass with respect to the standard basis. that are optimal in a least squares sense: they capture more power per mode than any other set of basis functions.

(We remark that this expansion is also known as the Karhunen-Loeve decomposition [START_REF] Tufillaro | Period doubling boundaries of a bouncing ball[END_REF].) The mean-square amplitude P, (or 'signal power') The optimality of the proper orthogonal modes therefore allows one to estimate the effective number of degrees of freedom excited in given motion by counting the minimum number of modes N such that

of
$1 3 CP,, (22)
for some fraction c of the total power. We now use the proper orthogonal decomposition to explore the relation between dimensionality, and spatial and temporal complexity in our impacting system. To apply the method, we take c = 0.99, and calculate the correlation matrix R, for various motions. The correlation matrix R,, was computed using 2 x 10" values of the configuration x(r) sampled at uniform intervals from the continuous-time response: the sampling rate chosen was in all cases roughly equal to the Nyquist frequency, although the results are not sensitive to this value. For each R,,, an error matrix E,, was also computed using the standard deviation from the mean for each component of R,,. The singular value decomposition was used to find the proper orthogonal modes 3, and the mean-square modal amplitudes A,. The largest singular value of E,, gives an upper bound on the error in the h,, which was used to ensure that sufficient data was taken so that the first N A, used in the calculations where above the computational noise floor.

The P,.,,, orbit of Figs 9 and 10 (o = 0.40816) furnishes an example of a nontrivial. yet spatially coherent, motion.

A portion of the evolution of the system in space and time is shown in Fig. 13. In Fig. 14 we show the first four proper orthogonal modes, and compare them to the normal modes for the system between impacts. The closeness of the proper orthogonal modes to the normal modes was typical of motions at all parameter values studied for this work, and suggests a simple interpretation of the cumulative statistical effect of repeated impacts: the spatial structure of the motion is close to what one would obtain by averaging over many impulse-response tests on the linear system. It must be emphasized, however, that the impacting system is strongly nonlinear (it is not linear even in the limit of 'small oscillations').

In Fig. 15, the cumulative power is plotted: one sees that only two proper orthogonal modes are needed to capture over 99% of the total power. We thus estimate that only 2 degrees of freedom are excited for this attractor. The correlation dimension dc was calculated using the Grassberger-Procaccia algorithm on a data set consisting of 10' points collected in the 20-dimensional impact Poincare section. For this and all other dimension estimates, the minimum absolute length scale considered valid was on the order of IO-', consistent with the cumulative accuracy of the time-step calculations: a characteristic 'knee' was observed on the correlation integral curves at this scale, indicating that scales below this level were in the system's noise floor. For the P,,,,, data, it was found that dc = 1.23 ? 0.01 over a scaling range of better than 1OOO:l. Thus. the P,,,,, orbit offers an example of a nonchaotic fractal attractor. Dimension theory [inequality (19)] predicts that the dimension m of the phase space is bounded by 2 and 5. which is consistent with the estimate of 2 proper orthogonal modes (2 phase-space dimensions per mode plus time yields a phase space with dimension 5). Figure 16(a) shows the spatial evolution for the chaotic orbit of Figs 11 and 12 (~0 = 0.35). Though somewhat more complex spatially than the orbit of Fig. 13, a high degree of spatial coherence is still visible. Indeed, one finds that only 3 proper orthogonal modes are required to exceed 99% of the total power. This is entirely consistent with the correlation dimension estimate of dC. = 2.21 f 0.04 over a 2O:l scaling range: dimension theory predicts a phase space dimension of between 3 and 7, and 3 modes plus time meet Arrows mark number at which cumulative power exceeds 99% of total (refer also to Table 1).

the upper bound exactly. A more strongly chaotic case is shown in Fig. 16(b), for which (r) = 0.36 and the damping parameter has been lowered to y = 0.01. The visually apparent spatial complexity is confirmed by the fact that 5 proper orthogonal modes are required to exceed 99% of the total power, and by the result dc = 6.2 k 0.2 (over a 2.2:1 scaling range). Five modes plus time give a phase space of dimension 11, which falls in the range predicted by dimension theory of 7 to 15. Table 1 displays additional information for these and other cases. The number of proper orthogonal modes N needed to exceed 99% of the total power is shown, along with the resulting estimate of the phase space dimension, 2 N + 1. These estimates should be compared to the bounds (19) obtained with the correlation dimension. Also included in the We remark that the resonance at w = 0.202 corresponds to a P,,, sticking motion, and the resonance at cc) = 0.408 corresponds to the P,,, sticking motion of Figs 9, 10, and 13. However, one should not conclude that all sticking periodic orbits lead to similar resonances, since several others exist at higher frequencies away from the resonances. Work aimed at a more precise understanding of this phenomenon is currently being carried out. Another feature of the bifurcation diagrams that requires further study is the banded pattern at higher frequencies in which a P,,, motion is alternated with a band of more complex responses (which include various P,,,, orbits as well as chaotic responses). Figure 6(b) shows that, for constant small input to the system, the energy level of the successive bands remains approximately constant.

Again, this pattern does not have an obvious relation to the natural frequencies of the non-impacting system. We have applied the proper orthogonal decomposition to study the relation between dimensionality and spatial structure for the system. A simple, total power criterion was used to estimate the number of active degrees of freedom for various motions (the results are summarized in Table 1). Clearly, taking the cutoff to be 99% of the total power is 19 somewhat arbitrary, and a more objective criterion would be preferable. Nevertheless, this simple method allowed us to examine the predictions of dimension theory in a consistent fashion: in five of the ten cases examined, the estimated number of phase-space dimensions (2N + 1) fell within the bounds predicted by inequalities (19). Four of these five were chaotic orbits (with (0 equal to 0.35, 0.36, 0.36, and 1.916), and even the remaining chaotic orbit (with w = 1.346) was reasonably close to the predicted bounds. The remaining result that fell within the predicted bounds was for a P,., sticking periodic orbit. The more regular orbits, however, failed to meet the predictions of dimension theory, with the poorest results belonging to the orbits which were simplest (in a topological sense): the P,,3 and PI., orbits for w equal to 0.870 and 1.750, respectively.

It is clear that rapid accelerations in such motions (which lead to the 'sharp edges' in the spatial evolution diagram of Fig. 17) necessitate the inclusion of more modes in any expansion: it may be that chaotic orbits tend to 'smooth out' such edges, thus explaining the better match with the predictions of dimension theory. At a more fundamental mechanical level, the singular nature of the impact condition tends to make the system dimensionally nonhomogeneous with respect to time: one should expect more modes to be excited 'just after' impact than on average. The uniform, continuous time average we did to obtain our proper orthogonal modes clearly has a big influence on the nature of the results. In contrast, preliminary work indicates that, as expected, averaging over impact Poincark section data, rather than continuous-time data, significantly increases the number of modes required. Current work is aimed at understanding the importance to the proper orthogonal decomposition method of such factors as sampling method, modal damping values. and impact-time statistics. We have found that the spatial structure, as evidenced by the proper orthogonal modes. is, in this case, virtually identical to that predicted by examination of the normal modes for the non-impacting system. This implies that the statistical effect of repeated impacts on the spatial structure of the nonlinear system is equivalent to what one would obtain by averaging over an ensemble of impulse response tests of the non-impacting linear system. Furthermore, this statistical equivalence does not appear to depend strongly on the type of response of the impacting system (in fact, we found much closer correspondence between the proper orthogonal modes and the normal modes for strongly chaotic responses). This suggests that deliberate excitation of chaotic impacting responses may be of practical use for engineering modal analysis. Of course. the proper orthogonal decomposition method is more general than standard modal analysis methods. since it makes no reference to linearity:

being entirely based on statistics, it will yield useful information even when the response between impacts is nonlinear.

Fig

  Fig. I. Configuration of the 10 degree of freedom impact oscillator

  1) u,(t;> = u,(t,,) = u,(t,,) (i = 1, 2,. . ., 9) (9.2) ~*o(CJ = ~(ur(t,J -uM(tJ) + UT(&) (9.3)

3 Fig. 2 .

 32 Fig. 2. Diagram showing the relation of the continuous-time system to the Poincare map defined by equations (10-12).

  (11.3) t ,,+ I = t,, + At,, (11.4)

  In Fig.3, we plot il, vs II during a sticking event. As shown in the figure, the A,* converge to the coefficient of restitution, right up to the limits of precision of the calculations. (A spurious divergence in Fig.3occurs after about 80 impacts. Assuming an initial time step of size O(l), this divergence occurs when At,, < (0.85)" = O(lO-'),

Fig. 3 .

 3 Fig. 3. Limiting behavior of impact time ratio during a sticking event. The ratios converge to the coefficient of restitution indicating that sticking occurs in finite time.

  Fig. 4. Bifurcation diagram showing relative tip velocity just after impact: computed with a fixed table amplitude A = 0. IO.

7 0Fig. 5 .Fig. 7 .

 757 Fig. 5. Blown-up portions of bifurcation diagram of Fig. 3. Three different sticking motions are labelled in the diagrams.

Fig. 8 .

 8 Fig. 8. Impact (top) and zero-phase stroboscopic (bottom) Poincark sections of the orbit shown in Fig. 7. In this and other figures, C#I = mod (ot, 2~) is the table phase.

Fig. 0 .

 0 Fig. 0. A P,,, orbit (co = 0.40816): (top) time series (x1,, is solid.table x , is dashed): (bottom) power spectrum.

Fig. 10 .

 10 Fig. 10. Impact (top) and zero-phase stroboscopic (bottom) Poincare sections of the orbit shown in Fig. Y.

Fig. 11 .Fig. 12 .

 1112 Fig. 11. A chaotic orbit (w = 0.35): (top) time series (xi,, is solid. table x , is dashed); (bottom) power spectrum.

  any motion is Ill P, = (x(r)*x(t)) = C(x,(r)q = &,(I,:, = ;A[.

Fig. 13 .Fig

 13 Fig. 13. Spatial evolution of the P,4 orbit of Figs 9 and 10 (2000 time steps with At = 0.2)

Fig. 16 .

 16 Fig. 16. Spatial evolution of chaotic orbits: (a) (I) = tl.35. y = 0.1; (b) w = 0.36, y = 0.01 (each has 2000 time steps with Ar = 0.2).

Fig. 17 .

 17 Fig. 17. Spatial evolution of the P1,3 orbit with w = 0.870 (2000 time steps with At = 0.1)

Table 1 .

 1 Summary of dimension and proper orthogonal decomposition results

		Y	Type of	N	2N + 1	4	[dcl	2(&j t I	y,,,*, 1 rrmn
			orbit
	0.34 0.35	0.1 0.1	p5.2 chaotic	3 3	7 7	1.00 f 0.00 2.21 2 0.04	1 3	3 7	220: I 20: 1
	0.36	0.1	h.32	3	7	1 .oo * 0.00	1	3	3.3:1
	0.36	0.01	chaotic	5	11	6.2 2 0.2	7	15	2.2:1
	0.36	0.001	chaotic	7	15	9.0 f 0.8	9-10	lY-21	1.8:1
	0.40816 0.870 1.336	0. 1 0.1 0.1	P = .i PI .? chaotic	2 5 3	5 11 9	1.23 IL 0.01 1.00* 2.58 f 0.05	2 1 3	5 3 7	1000: I _* 6.7: I
	1.750	0.1	PI h	5	11	1 .oo*	1	3	_*
	1.916	0. 1	chaotic	3	7	2.9 t 0.1	3	7	2.7:1

*No error estimates or scaling regions are presented for these frequencies since the impact Poincare section contained only one cluster of points, and hence all length scales were in the noise floor.

table is

 is especially for the chaotic cases, but fails, sometimes in spectacular fashion, for some of the periodic motions. Some insight into how this can happen can be found by examining Fig.17, in which the spatial evolution of the P1,3 response with w = 0.870 is shown. The sharp edges visible in the diagram,

	the size of the scaling region associated	with each dimension	estimate.	The table
	clearly shows that dimension	theory as summarized	by inequality	(19) performs reasonably
	well, which indicate	that the masses are being rapidly
	accelerated,	make it necessary to include more proper orthogonal	modes than required on
	the basis of purely topological	arguments.
						6. SUMMARY AND DISCUSSION
	In this paper we have presented	a numerical	study of bifurcation,	dimensionality,	and
	spatial coherence	in a 10 degree of freedom impact oscillator.	Bifurcation	diagrams for the
	system were obtained	over the entire range of natural frequencies	of the non-impacting
	system.	Several	characteristics	of the diagrams	merit further	investigation.	We have
	demonstrated,
							The energy bifurcation	diagrams	of Fig. 6 show two
	nonlinear	resonances	which do not occur at non-impacting	natural frequencies:	nor are the
	resonances		related		to an obvious,	simple combination	resonance.

via relations (17) and (18), and in Fig.

6(b)

, that the input to the system is controlled by the table velocity.

degree of freedom impact oscillator

Ach-nolc'le~Re,,zerIts-This work wa5 supported by the National Science Foundation grant MSS 8Y-10814. and by the Air Force Office of Scientific Reacarch grant AFOSR-Yl-0330.

Numerical

Recipes: The Art of Scientific Computing.

Cambridge

University Press, New York (1986)

Stochastic Finite Elements:

A Spectral Approach.

Springer, New York (1991). APPENDIX When the system of Fig. 1 is between impacts, it is governed by equations (I) and ( 2). The natural frequencies of the non-impacting system are found by standard numerical means to be as shown in Table Al.


D J Wagg

Periodic sticking motion in a two-degree-of-freedom impact oscillator

Keywords: Impact, 2DOF Oscillator, Periodic, Sticking, Multi-sliding

Periodic sticking motions can occur in vibro-impact systems for certain parameter ranges. When the coefficient of restitution is low (or zero), the range of periodic sticking motions can become large. In this work the dynamics of periodic sticking orbits with both zero and non-zero coefficient of restitution are considered. The dynamics of the periodic orbit is simulated as the forcing frequency of the system is varied. In particular, the loci of Poincaré fixed points in the sticking plane are computed as the forcing frequency of the system is varied. For zero coefficient of restitution, the size of the sticking region for a particular choice of parameters appears to be maximized. We consider this idea by computing the sticking region for zero and nonzero coefficient of restitution values. It has been shown that periodic sticking orbits can bifurcate via the rising/multi-sliding bifurcation. In the final part of this paper, we describe three types of post-bifurcation behavior which occur for the zero coefficient of restitution case. This includes two types of rising bifurcation and a border orbit crossing event.

Introduction

In this paper periodic sticking motions which occur in the dynamics of a two-degree-of-freedom impact oscillator are considered. The impact oscillator consists of two masses, coupled with springs and dashpots, and the motion of both the masses is restricted by rigid constraints [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF]. These type of systems can be used to model a range of physical applications mainly in mechanical engineering [START_REF] Masri | Theory of the dynamic vibration neutraliser with motion-limiting stops[END_REF][START_REF] Moon | Chaotic vibrations of a beam with non-linear boundary conditions[END_REF][START_REF] Chatterjee | On impact dampers for non-linear vibrating systems[END_REF][START_REF] Bapat | The general motion of an inclined impact damper with friction[END_REF][START_REF] Blazejczyk-Okolewska | Analysis of an impact damper of vibrations[END_REF][START_REF] Semercigil | Experiments with tuned absorber-impact damper combination[END_REF][START_REF] Azeez | Numerical and experimental analysis of a continuous overhung rotor undergoing vibroimpacts[END_REF][START_REF] Metallidis | Vibration of a continuous system with clearance and motion constraints[END_REF]. Many mathematical-numerical studies of these systems have been carried out, and particular interest has been focused on bifurcation behavior [START_REF] Shaw | The onset of chaos in a two-degree of freedom impacting system[END_REF][START_REF] Luo | Hopf bifurcation of a two-degreeof-freedom vibro-impact system[END_REF][START_REF] Wen | Codimension-2 hopf bifurcation of a two-degreeof-freedom vibro-impact system[END_REF][START_REF] Luo | Hopf bifurcations and chaos of a two-degree-of-freedom vibro-impact system in two strong resonance cases[END_REF][START_REF] Hogan | Graph theory and piecewise smooth dynamical systems of arbitary dimension[END_REF][START_REF] Lenci | Regular nonlinear dynamics and bifurcations of an impacting system under general periodic excitation[END_REF][START_REF] Valente | Hybrid dynamics of two coupled oscillators that can impact a fixed stop[END_REF]. More general studies of multiple degree of freedom impact oscillators have also been carried out [START_REF] Theodossiades | Periodic and chaotic dynamics of motor-driven gear-pair system with backlash[END_REF][START_REF] Cusumano | Period-infinity periodic motions, chaos and spatial coherence in a 10 degree of freedom impact oscillator[END_REF][START_REF] Nigm | Effect of an impact damper on a multi-degree of freedom system[END_REF][START_REF] Pfeiffer | Multibody Dynamics with Unilateral Contacts[END_REF][START_REF] Natsiavas | Dynamics of multiple-degree-of-freedom oscillators with colliding components[END_REF][START_REF] Pun | Forced vibration of a multidegree impact oscillator[END_REF], but these focus primarily on a single impact constraint.

In mechanical systems with vibration and impact, chatter and sticking are phenomena which have been observed for a wide range of parameter values. Chatter and sticking in single degree of freedom impact oscillators has been studied in detail [START_REF] Budd | Chattering and related behaviour in impact oscillators[END_REF][START_REF] Nordmark | On chattering bifurcations in 1dof impact oscillator models[END_REF] and also noted to occur in two-degree-of-freedom systems, particularly for low forcing frequencies [START_REF] Wagg | Chatter, sticking and chaotic impacting motion in a two-degree of freedom impact oscillator[END_REF]. Vibro-impact systems exhibit a rich variety of periodic motions (see for example [START_REF] Bapat | Periodic motions of an impact oscillator[END_REF]), and periodic sticking motions can be found for particular parameter values in both single and multi-degree-of-freedom systems. Multi-degreeof-freedom systems with a single constraint have been studied by Toulemonde and Gontier [START_REF] Toulemonde | Sticking motions of impact oscillators[END_REF] and Wagg and Bishop [START_REF] Wagg | Chatter, sticking and chaotic impacting motion in a two-degree of freedom impact oscillator[END_REF], where periodic sticking motions were observed for low forcing frequencies. In [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF], periodic sticking motions were noted for a two-degree-offreedom impact oscillator with two constraints, and it was shown how the entry boundary to the sticking region can be defined. The rising bifurcation discussed by Toulemonde and Gontier [START_REF] Toulemonde | Sticking motions of impact oscillators[END_REF] has been shown to be equivalent to the multi-sliding bifurcation [START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF], which occurs in the study of nonsmooth systems [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF][START_REF] Kowalczyk | On a novel class of bifurcations in hybrid dynamical systems[END_REF][START_REF] Bernardo | Bifurcations of dynamical systems with sliding: derivation of normal-form mappings[END_REF].

In this paper, we study the dynamics of periodic sticking motions which exist in a two-degree-offreedom impact oscillator with motion constraints on both masses [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF]. For low (or zero) coefficient of restitution, the range of periodic sticking motions becomes large. This is significant because the maximum extent of the sticking region should be defined for the zero case. This idea is considered by using a comparison of the sticking region obtained by simulating the dynamics of periodic sticking orbits with both zero and non-zero coefficient of restitution. Then the loci of Poincaré fixed points in the sticking plane are computed as the forcing frequency of the system is varied. Projections of these loci into the system state space indicate that the zero case does not define the largest region of sticking motions for the examples considered.

In the final part of this paper, we consider the postbifurcation behavior following rising bifurcations in the zero coefficient of restitution case. We first illustrate the 'standard' rising bifurcation behavior described by Wagg and Bishop [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF][START_REF] Toulemonde | Sticking motions of impact oscillators[END_REF] and following this we describe one other example where the rising has a receding behavior. Finally we discuss an example were a sticking orbit passes through a border orbit [START_REF] Nordmark | On chattering bifurcations in 1dof impact oscillator models[END_REF], which defines the limiting extent of the sticking region.

Mathematical model

We consider a coupled two-degree-of-freedom system, which is shown schematically in Fig. 1. This sys-

k 1 k 2 s 2 f 1 f 2 x 1 x 2 s 1 c 1 c 2 m 1 m 2
Fig. 1. Schematic representation of an N degree of freedom impact oscillator with multiple motion limiting constraints.

tem has already been described in detail in [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF][START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF]. The governing equations for the system away from impact can be expressed as

m 1 ẍ1 + c 1 ẋ1 + c 2 ( ẋ1 -ẋ2 ) + k 1 x 1 + k 2 (x 1 -x 2 ) = f 1 , (1)
m 2 ẍ2 + c 2 ( ẋ2 -ẋ1 ) + k 2 (x 2 -x 1 ) = f 2 , (2)
where x 1 represents the displacement of mass m 1 and x 2 the displacement of mass m 2 . The spring stiffnesses are given by k 1 , k 2 and the damping constants by c 1 , c 2 and the distance to the motion constraints are given by s 1 and s 2 , respectively. The harmonic forcing functions are f 1 = A 1 cos( t) and f 2 = A 2 cos( t). Eq. ( 2) has a dual condition for free flight that (x i -s i ) < 0 for s i > 0 and (x i -s i ) > 0 for s i < 0 which can be written as (x i -s i )≶0 ∀s i ≷0 for i = 1, 2. We also assume that the distance between masses is large enough so that they do not impact with each other. Eqs. ( 1) and ( 2) can be written in the nondimensionalized form

m 0 0 1 ¨ 1 ¨ 2 + 2 1 √ m k + 2 2 -2 2 -2 2 2 2 ˙ 1 ˙ 2 + 1 + k -1 -1 1 1 2 = f1 f2 , ( 3
)
where m = m 1 /m 2 , k = k 1 /k 2 , 1 = c 1 /(2m 1 n1 ), 2 = c 2 /(2m 2 n2 ), n1 = √ k 1 /m 1 , n2 = √ k 2 /m 2 , = / n2 , f1 = P 1 cos( ), f1 = P 2 cos( ), P 1 = A 1 /(k 2 x c ), P 2 = A 2 /(k 2 x c ), = n2 t and = x/x c .
The non-dimensional variable is achieved by dividing displacement, x, by a constant displacement x c . This choice is arbitrary, and therefore we will assume that x c = 1, such that the non-dimensional distances to the motion constraints are 1 = s 1 /x c and 2 = s 2 /x c . The non-dimensional phase, , is defined as = mod2 / .

The parameter values have been selected as m 1 = m 2 = k 1 = k 2 = 1 and c 1 = c 2 = 0.1 which means that in the non-dimensionalized case m = k = n1 = n2 = 1 and 1 = 2 = . These specific parameters were chosen to give a simple relationship between the natural frequencies and the system eigenvalues, and the damping value is chosen to represent a physically realistic choice for a mechanical spring-mass-damper system. In this case, Eq. (3) simplifies to

¨ 1 ¨ 2 + 2 2 -1 -1 1 ˙ 1 ˙ 2 + 2 -1 -1 1 1 2 = P 1 cos( ) P 2 cos( ) . ( 4
)
The natural frequencies of the non-dimensional system are given by j for j = 1, 2 where j are the eigenvalues of the 2 × 2 coupling matrix

[E] = 2 -1 -1 1 . (5)
The corresponding normalized eigenvectors j can be used to construct a orthogonal modal matrix

[ ] = [{ 1 }, { 2 }].
We can then transform Eq. ( 4) into a modal form by defining modal coordinates q={q 1 , q 2 } T , such that = [ ]q and

[I ] q + 2 q + [ ]q = [ ] T f (t), (6)
where

[ ] = [ ] T [E]
[ ] is the diagonal matrix of the eigenvalues, j , j = 1, 2 and f (t) = [P 1 cos( ), P 2 cos( )] T . In this modal formulation, we define the vector i = [ i1 , i2 ], such that an impact occurs when i q = i , i = 1, 2. Hence Eq. ( 6) is valid only for ( i qi )≶0 ∀ i ≷0, which is equivalent to the condition that (x i -s i )≶0 ∀s i ≷0 for the ith impacting mass. For this system there are two modal impact vectors, 1 = [ 11 , 12 ] and 2 = [ 21 , 22 ], such that at impact 1 q = 1 and 2 q = 2 , where q = [q 1 , q 2 ] T .

Eq. ( 6) for each mode (with P 2 = 0) is given by qj + 2 j j qj + j q j = 1i P 1 cos( ), j = 1, 2.

Eq. ( 7) has a well-known exact solution for underdamped oscillations 0 < j < 1 [START_REF] Timoshenko | Vibration Problems in Engineering[END_REF], such that for each mode exact solutions can be found between impacts [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF]. We consider only the under-damped case as this is the case of most interest for mechanical systems.

For the numerical simulations in this paper we set the forcing amplitudes as P 2 = 0 and P 1 = 0.5 and take initial conditions q 1 (t 0 ) = q 2 (t 0 ) = q1 (t 0 ) = q1 (t 0 ) = t 0 = 0. When ( ii ) = 0 for i = 1, 2 an impact occurs which is modelled using an instantaneous coefficient of restitution rule [START_REF] Thompson | Nonlinear Dynamics and Chaos[END_REF] . For single impacts the coefficient of restitution rule is ẋi (t + ) = -r ẋi (t -), where, t -is the time just before impact, t + is the time just after impact and r is the coefficient of restitution with a value in the range r ∈ [0, 1]. In matrix form the coefficient of restitution rule can be written as

˙ (t + ) = [R k ] ˙ (t -)
where for the system being considered there are three different cases for the

[R k ] matrices [R 1 ] = -r 0 0 1 , [R 2 ] = 1 0 0 -r , [R 3 ] = -r 0 0 -r , (8)
corresponding to mass 1 impacting, mass 2 impacting and simultaneous impact of both masses.

In modal form the coefficient of restitution rule becomes [ ] q( + ) = [R k ][ ] q( -), which leads to the relation for the modal velocities after impact of q( + ) = [ Rk ] q( -), where

[ Rk ] = [ ] -1 [R k ][ ]
is the set of matrices which represents a linear transform of modal velocities just before impact to modal velocities just after impact for the three possible impact cases [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF].

Sticking solutions

For this system there are two possible sticking regimes; when 1 = 1 and when 2 = 2 . Each regime has a reduced set of governing equations with explicit solutions [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF]. It is also possible to have a dual sticking regime when both 1 = 1 and 2 = 2 simultaneously, with the result that there are no dynamics in the system.

In the case where mass 1 sticks 1 = 1 and ˙ 1 = 0, and the equations of motion reduce to a single equation

¨ 2 + 2 ˙ 2 + 2 = 1 . ( 9
)
The condition for the mass being held in place during sticking is related to the force on the mass. The equivalent non-dimensional release expression is given by

F 2 = 2 ˙ 2 + 2 -2 1 + P 1 cos( ). ( 10
)
The forcing term P 2 = 0 in Eq. ( 9), but P 1 occurs in the release expression, so this case is referred to as the forced sticking case. The sticking phase ends when F 2 becomes zero and changes sign at which time = f . Eq. ( 9) has the exact solution

2 = e -( -s ) (C 1 cos[ 1 -2 ( -s )] + C 2 sin[ 1 -2 ( -s )]) + 1 . ( 11
)
At the start of the sticking period s = and the constants C 1 and C 2 are found using the initial conditions for 2 and ˙ 2 when 1 ( s ) = 1 and

˙ 1 ( s ) = 0 [1].
In the case 2 = 2 and ˙ 2 = 0, the reduced equation of motion is given by

¨ 1 + 2 ˙ 1 + 2 1 -2 = P 1 cos( t). ( 12
)
The release condition is governed by

F 1 = ˙ 1 + 1 -2 . ( 13
)
This is therefore referred to as the unforced sticking case.

Eq. ( 12) has the exact solution

1 = e -2 ( -s ) (C 1 cos(2 1 -2 ( -s )) + C 2 sin(2 1 -2 ( -s ))) + C 3 cos( t -) + 2 /2, ( 14
)
where = arctan((4

)(2 -2 )).
The initial conditions for both reduced equations can be taken directly from the appropriate values of i and ˙ i immediately prior to a sticking phase. These initial conditions allow the computation the constants C 1 , C 2 prior to the beginning of the next phase of motion, and C 3 is found as part of the particular solution [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF].

Periodic sticking motion

An example of the type of periodic sticking orbit which will be considered in this paper is shown schematically in Fig. sticking occurs for mass 1 between A and A , and for mass 2 between B and B . One complete period occurs between points A and C. From A to A mass 1 is stuck and mass 2 is in free flight. Then from A to B both masses are in free flight, followed by a second sticking regime from B to B where mass 2 is stuck and mass 1 is in free flight. Finally, both masses are in free flight between B and C. The period of (nondimensional) time spent sticking is T si for i = 1, 2 and for periodic sticking orbits the proportion of the whole period spent sticking is p i = T si /T , where T = 2 / so that p i = T si /2 .

The change from free motion of both masses to one mass sticking represents a reduction in the degree of freedom of the system from 2 to 1. For the example in Fig. 2 this reduction from 2 to 1 occurs at points A and B. Similarly, there is an increase in degree of freedom from 1 to 2 at points A and B .

The two-degree-of-freedom system has a parameter set = { , P 1 , P 2 , , 1 , 2 , r}. For any particular choice of these parameter values there could typically be regions of non-impacting behavior, vibro-impacting behavior, chatter and sticking [START_REF] Wagg | Chatter, sticking and chaotic impacting motion in a two-degree of freedom impact oscillator[END_REF]. We will assume that P 2 =0, and that , P 1 , , 1 , 2 have suitable constant values such that for a range of excitation frequencies , vibro-impacting motions, chatter and sticking occur. For such a choice of parameters, the coefficient of restitution, r, will then define the extent of the potential sticking region, S, for periodic sticking orbits.

For r = 0, the impacts are completely plastic, and sticking will generally occur after every impact-in Section 4.1 we show an exception to this where impact occurs without sticking. For r = 0, sticking will only occur after a complete chatter sequence [START_REF] Budd | Chattering and related behaviour in impact oscillators[END_REF], and as a result sticking orbits will occur for lower excitation frequencies-because this allows chatter to become complete. As r → 1 impacts tend towards being perfectly elastic and sticking orbits will not appear even for low excitation frequencies. A numerical example is shown in Fig. 3, where we compare two time series. In Fig. 3(a), the coefficient of restitution is r =0.7, and in Fig. 3(b), the coefficient of restitution is r = 0.0. It is clear that for the r = 0.0 the sticking proportions of the periodic orbit, p 1 , p 2 are at their maximum values for all the other parameters in fixed.

Classifying periodic sticking orbits

Periodic impacting orbits are usually classified by the number of forcing periods, n, and the number of impacts, m, which occur during one period of the motion, denoted P (n, m) [START_REF] Shaw | A periodically forced piecewise linear oscillator[END_REF]. Periodic orbits with sticking have an infinite number of instantaneous impacts during one period of motion [START_REF] Budd | Chattering and related behaviour in impact oscillators[END_REF], so every periodic sticking orbit is classified P (n, ∞). With the understanding of multi-sliding bifurcations [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF][START_REF] Kowalczyk | On a novel class of bifurcations in hybrid dynamical systems[END_REF][START_REF] Bernardo | Bifurcations of dynamical systems with sliding: derivation of normal-form mappings[END_REF], and how they occur in impacting systems [START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF], there is a further distinction which can be made for periodic sticking orbits-the number (and duration) of the sticking phase(s). For example, if k is the number of sticking phases per period we could classify the orbit using P (n, ∞, k). A more detailed classification would include the proportion of the period spent sticking, p =[p 1 , p 2 ], such that P (n, ∞, k, p). In fact the majority of periodic sticking orbits considered in this work are P (1, ∞, 1). In the case where a periodic motion occurs with a sticking phase and an additional number of separate impacts we use the notation P (n, ∞ + m, k), an example of this is discussed in Section 4.1.

The sticking region

Sticking behavior in vibro-impact systems is analogous to sliding behavior in some electrical systems [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF], as discussed in [START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF]. For studying the behavior of sticking orbits in general, it is useful to define the region in the system phase space where these orbits exist. The sticking/sliding region(s) consists of a manifold(s) within the system state space on which the sticking/sliding orbits exist. These sticking/sliding orbits are then (usually) restricted to some region on the manifold by conditions which define the entry (start of sliding) and exit (end of sliding) of the orbit from the manifold. For some systems [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF] it is possible to define the manifold and the entry and exit boundaries explicitly. For vibro-impact systems is has been pointed out that only the manifold and exit boundary can be defined explicitly [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF].

For the system studied here there are two manifolds in the system state space on which sticking can take place. These are defined by the impact conditions of the system such that if = { 1 , 2 , ˙ 1 , ˙ 2 , } T is the state vector, then the system state space can be defined as G = { ∈ R 4 × : 1 1 , 2 2 }. Then the two impact manifolds are defined as 1 = { ∈ R 4 × : sticking orbit to leave the sticking region S i is given by the equations for F i = 0, for i = 1, 2-Eqs. [START_REF] Shaw | The onset of chaos in a two-degree of freedom impacting system[END_REF] and [START_REF] Luo | Hopf bifurcations and chaos of a two-degree-of-freedom vibro-impact system in two strong resonance cases[END_REF].Therefore, we can analytically define the exit boundary of each sticking region jS i , using the conditions for F i = 0, for i = 1, 2 in Eqs. ( 10) and [START_REF] Luo | Hopf bifurcations and chaos of a two-degree-of-freedom vibro-impact system in two strong resonance cases[END_REF].

Although there is no way of analytically defining the entry boundary of the sticking regions, it would appear from Fig. 3 that the entry boundary could be defined by the point of first sticking when r =0, i.e. the maximum possible extent of the sticking region occurs when no chatter precedes the sticking. Therefore, to test whether the r = 0 case does define the maximum possible extent of the sticking region, the loci of the first sticking points for mass 2 sticking (point B in Fig. 2) have been computed as the forcing frequency is varied. Because mass 2 is sticking, the state variables which will define the sticking orbit are 1 , ˙ 1 and , and we will plot these values at the point where mass 2 first sticks, i.e. the entry point into the sticking region. For the r = 0 case the bifurcation diagrams showing the amplitude of 1 and at the first sticking point as is varied through the range 0.1-2.5 are shown in Figs. 4(a) and(c), respectively. As a comparison we have shown the same plots for the r =0.7 case in Figs. 4(b) and(d), from which we note that sticking exists for a much larger range of values in the r = 0 case. Note also that there are some non-smooth jumps in the variation of with shown in Figs. 4(c) and(d). This is explained by observing the modulo value, 2 / , which is plotted as a dashed line in Figs. 4(c) and(d). For both cases shown in Fig. 4, r = 0 (Figs. 4(c)) and r = 0.7 (Fig. 4(d)), the loci begins at a value below the modulo value. However in both cases, as is increased, the loci eventually intersect with the modulo line such that the phase value becomes zero. This explains the jumps of the loci in Figs. 4(c)and(d), at ≈ 0.75 for Fig. 4(c) and ≈ 0.21 in Fig. 4(d). This also accounts for the division of the loci into two sections in Figs. 5(c) and(d).

As is varied the first sticking points defined by 1 , ˙ 1 and describe a locus in R 3 . In Fig. 5(a) the projection of this locus into the 1 , ˙ 1 plane is shown, and in Fig. 5(c) the projection into the , ˙ 1 plane is shown. These loci define the entry into the sticking region, and the exit boundary jS 1 is marked as a dashed line in Fig. 5(a).

By comparing the zero and non-zero case, we can see that the loci have a particularly complex structure when projected into the 1 , ˙ 1 plane, Figs. 5(a) and(b). It is also clear that the r = 0 loci does not encompass the entire region occupied by the r = 0.7 loci for the 1 , ˙ 1 projection. However, the explanation of this is apparent from the time series in Fig. 3, where it can be seen that due to the time at which the sticking starts, the magnitude of 1 value is greater for the r = 0.7 case than for the r = 0 case. This can also be seen by comparing Figs. 4(a) and(b) where the value of 1 at the first sticking point is plotted against . For the , ˙ 1 projections shown in Figs. 5(c) and (d), it can be seen that the r = 0 loci does not encompass the entire region occupied by the r = 0.7 loci, but it does have higher values-indicating a stronger possibility for defining a bound on the region S 1 .

It is interesting to note that as the system parameter is varied smoothly, the resulting first sticking loci contain several discontinuities-non-smooth points and discontinuous jumps. Apart from phase transitions, these points represent non-smooth changes (bifurcations) in the dynamics of the periodic orbits, and will be discussed in Section 4.

Sliding bifurcations

The sliding orbits in electrical systems have been shown to exhibit particular types of sliding bifurcations under parameter variation [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF]. There are four types of sliding bifurcation which can occur [START_REF] Kowalczyk | On a novel class of bifurcations in hybrid dynamical systems[END_REF], for which the normal form mappings have been derived [START_REF] Bernardo | Bifurcations of dynamical systems with sliding: derivation of normal-form mappings[END_REF]. A multi-sliding bifurcation is one of the four cases, which occurs in the systems studied by Di Bernardo et al. [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF], Kowalczyk and di Bernardo [START_REF] Kowalczyk | On a novel class of bifurcations in hybrid dynamical systems[END_REF], and Di Bernardo et al. [START_REF] Bernardo | Bifurcations of dynamical systems with sliding: derivation of normal-form mappings[END_REF], and is the most significant for our current study. Previous physical examples of the multi-sliding bifurcation have been studied in models of relay feedback systems [START_REF] Bernardo | Self-oscillations and sliding in relay feedback systems: symmetry and bifurcations[END_REF] and friction oscillators [START_REF] Bernardo | Sliding bifurcations: a novel mechanism for the onset of chaos in dry friction oscillators[END_REF].

In addition to multi-sliding, in the case of r = 0, the system will normally have a grazing-sliding bifurcation (in this example grazing followed by sticking) each time a grazing event occurs. This is similar to the grazing-sliding in the friction oscillator example studied by Di Bernardo et al. [START_REF] Bernardo | Sliding bifurcations: a novel mechanism for the onset of chaos in dry friction oscillators[END_REF], but will not be considered in detail here.

Three examples of post-bifurcation behavior

The multi-sliding (or rising) bifurcation occurs when a sliding (sticking) orbit touches the boundary of the sliding (sticking) region S i [START_REF] Toulemonde | Sticking motions of impact oscillators[END_REF][START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF]. Physically this means that the force holding the mass against the constraint becomes zero, and as the bifurcation parameter continues to vary, the mass lifts off (or rises) from the constraint. This results in a sudden reduction in sticking time [START_REF] Toulemonde | Sticking motions of impact oscillators[END_REF][START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF].

In this subsection, we look at three different types of post-bifurcation behavior which relate to observations of the behavior of the first sticking loci shown in Figs. 4 and5 for the r = 0 case. All three examples are for mass 1 sticking. The first example is shown in Fig. 6, and occurs close to = 0.254. In Fig. 6(a) and (b) we see that a rising has occurred some where near the middle of the sticking phase. As increases the rise propagates towards the release point for the mass, Fig. 6(c) and (d). The result is that the single sticking region is divided into two parts. So the periodic sticking orbit goes from P (1, ∞, 1) → P (1, ∞, 2). However the sticking phase following the rise quickly decreases until a limit point, where it becomes a single point (like grazing) just before the mass lifts off-Fig. 6(d) then (f). This type of post multi-sliding behavior will be called the standard case, and has been discussed by Toulemande and Gontier [START_REF] Toulemonde | Sticking motions of impact oscillators[END_REF], Wagg [START_REF] Wagg | Rising phenomena and the multi-sliding bifurcation in a two-degree of freedom impact oscillator[END_REF].

A second example of post multi-sliding behavior is shown in Fig. 7, and occurs close to = 0.475. In this case the rise is very close to the first sticking point. As with the standard case, the periodic sticking orbit goes from P (1, ∞, 1) → P (1, ∞, 2), but one of the new sticking phases is very small. The post-multisliding behavior is then that the rise becomes larger in amplitude, before reaching a maximum and then declining until P (1, ∞, 2) → P (1, ∞, 1) through a reverse multi-sliding event. Note, this event can be seen in Fig. 4(a) as a series of points below the main loci close to = 0.475. This is called the receding multi-sliding case.

The last example, shown in Fig. 8, occurs for frequency values close to = 0.7. In this example the sliding orbit passes through what has now been termed a border orbit-the orbit on which the sticking zone shrinks to zero [START_REF] Nordmark | On chattering bifurcations in 1dof impact oscillator models[END_REF]. In this case, the mass first impacts without sticking, and then impacts with a sticking phase, classified as a P (1, 1 + ∞, 1) periodic orbit. It is worth noting that in this type of multi-mass system it is possible for an impact to occur without sticking even for the r = 0 case. This can be seen from Eq. [START_REF] Shaw | The onset of chaos in a two-degree of freedom impacting system[END_REF], which is the release condition for when mass 1 is sticking. When an impact occurs, Eq. ( 10) must have the same sign as the constraint distance (negative in this case) for sticking to occur. In effect, Eq. ( 10) is the non-dimensional equivalent of the force holding the mass against the constraint [START_REF] Wagg | Dynamics of a two degree of freedom vibro-impact system with multiple motion limiting constraints[END_REF]-if the mass impacts and the force is acting away from the constraint, then no sticking will occur. It is clear that this force is dependent on the displacement and velocity of mass two and the external forcing. For example, computing the values of F 2 for Fig. 8(b) ( = 0.85), at the first impact ( =964.12) F 2 =0.202, so no sticking occurs. At the second impact ( = 965.65) F 2 = -0.1594, so sticking does occurs.

As is increased, the sticking phase (for mass 1) reduces to zero, until the border orbit is reached (Fig. 8(c)). Beyond this, mass 1 has no sticking phase (although mass 2 continues to stick until a border orbit close to = 1.4-Fig. 4(a)). The border orbit event for mass 1 has a clear effect on the stick-ing values for mass 2 which can be seen clearly in Fig. 4(a) as the discontinuity close to = 0.9. This effect can also be seen clearly as a sharp discontinuity in Figs. crossing the periodic orbit makes the transition from P (1, 1 + ∞, 1) → P (1, 1)-sticking periodic orbit to impacting periodic orbit.

Conclusions

In this paper, we have examined the behavior of periodic sticking orbits which occur in a two-degreeof-freedom impact oscillator. For these periodic sticking orbits, two cases have been considered, one with r = 0.7 and the other with r = 0. For each case we have computed the loci of first sticking points in the sticking region-as forcing frequency is varied-demonstrating the complex nature of the entry boundary for this region. The projections of the sticking region into the 1 , ˙ 1 , and , ˙ 1 planes indicated that the condition r = 0 did not bound the stick-ing region completely-although in the , ˙ 1 plane the first sticking loci provided a partial boundary to the region.

Plots of the loci of first sticking points against the bifurcation parameter (forcing frequency) showed clear non-smooth jumps and discontinuities. Three of these non-smooth events have been discussed in more detail, including two types of multi-sliding bifurcation and a border orbit crossing. The discovery of multi-sliding bifurcations and border orbit crossing in models of mechanical and electrical systems is a very recent addition to the literature on this subject. In this paper we have shown, by example, how border orbit crossing events manifest themselves. We have also shown that multi-sliding bifurcations can have two distinct types of post bifurcation behavior-the standard case and the receding case involving a reverse multi-sliding event.

Fig. 2 .

 2 Fig. 2. Schematic diagram of a periodic sticking orbit. Sticking occurs between A to A and B to B . Solid line mass 1, dashed line mass 2.

Fig. 3 .

 3 Fig. 3. Time series showing; (a) r = 0.7 and (b) r = 0. The maximum potential sticking region occurs when r = 0. Parameters = 0.05, P 1 = 0.5, = 0.4, 1 = -0.3, 2 = 0.1.

Fig. 4 .

 4 Fig. 4. First sticking loci as is varied: (a) and (b) values of 1 at the point when mass 2 sticks to the constraint. (c) and (d) corresponding values of . (a) and (c) r = 0, (b) and (d) r = 0.7. Parameters = 0.05, P 1 = 0.5, 1 = -0.3, 2 = 0.1.

Fig. 5 .

 5 Fig. 5. First sticking loci as is varied from 0.1 to 2.5: (a) and (b) projected onto the ˙ , plane; (c) and (d) projection onto the ˙ , plane; (a) and (c) r = 0.0, (b) and (d) r = 0.7. The dashed line in (a) and (b) denotes the boundary of the sticking region jS where the constraining force drops to zero, and the mass is released from sticking. Parameters = 0.05, P 1 = 0.5, 1 = -0.3, 2 = 0.1.

Fig. 6 .

 6 Fig. 6. Rising bifurcation sequence 1: Solid line is mass 1, dashed line mass 2. Parameters = 0.05, P 1 = 0.5, 1 = -0.3, 2 = 0.1, r = 0; (a) and (b) = 0.254, (c) and (d) = 0.26, (e) and (f) = 0.265.

Fig. 7 .

 7 Fig. 7. Rising bifurcation sequence 2: Solid line is mass 1, crosses are computation points. Parameters =0.05, P 1 =0.5, 1 =-0.3, 2 =0.1, r = 0, = 0.469-0.505.

Fig. 8 .

 8 Fig. 8. Border orbit crossing: Solid line is mass 1. Parameters = 0.05, P 1 = 0.5, 1 = -0.3, 2 = 0.1, r = 0, = 0.7 -0.91.

= 1 , ˙ 1 =0} and

={ ∈ R 4 × : 2 = 2 , ˙ 2 =0}.On each impact manifold, i , a corresponding sticking region S i exists for i = 1, 2. The condition for a

Acknowledgements

This work was supported as part of an Advanced Research Fellowship from the EPSRC. The author would like to thank Piotr Kowalczyk for constructive comments regarding this work.