Stable periodic vibroimpacts of an oscillator

Abstract The asymptotically stable vibrations of a loaded oscillator colliding periodically with a rigid mass are described. Comparison of the numerical results with the few existing examples is encouraging but inconclusive. Better overall agreement is demonstrated with fairly comprehensive measurements from a specially built experimental rig. Impact motions are shown to be very sensitive to small fluctuations in the clearance between masses and the stiffness and loading of the oscillator near its linear, or collisionless, resonant frequency. The rigid mass is quite an effective damper at or just above this frequency condition.


INTRODUCTION
A general theory to describe the stable periodic motions of two colliding rigid bodies called an impact-pair has been developed and verified experimentally in a companion paper [l]. The theory will be extended by letting one of the bodies have flexibility and damping in order to assess resonant effects. Such effects are likely important in high speed mechanisms or machinery where clearances arise from tolerances, wear, etc. On the other hand, an extra impacting body may be introduced deliberately into a resonant system to reduce excessive vibration excursions. The device is normally termed an "impact damper" and can be quite effective despite the high transient accelerations generated by the impacts [2-211. The basic theory of these applications will be checked against previous, rather restricted, theoretical results and new, exhaustive experimental data. A closed form solution has been derived for impacts consistently happening twice for every cycle of a sinusoidal load [2-5, 7, 8, 10, 11, 14, 22, 231. The impacts have been assumed generally instantaneous and distributed primarily evenly in time. Masri [7, 131 and the Kobrinskii's [23] subsequently applied the concept of error propagation in difference equations to ascertain whether these periodic motions were asymptotically stable or not. These approaches are fundamental but the practical application of the solutions is quite restricted. For example Sadek [6] and Dittrich [9] suggested from limited experimental and theoretical studies that uneven temporal distributions are most likely at frequencies close to where the collisionless, flexible system is resonant.
Sadek also inferred after careful experimentation that "equally spaced impacts hardly ever occur" even away from this resonance.
Therefore it is not surprising that Sadek [ 1 cycle were emphasized; only a single, probably analog computed illustration of three impacts was given by Masri alone. The lone three impact example may well be the consequence of a slowly converging and, hence, computationally time consuming iteration scheme for non-linear equations. An efficient algorithm will be used here for a comparable set of linearized equations to facilitate the computation of stable periodic motions involving any number of impacts. Numerical examples however will be confined to five or fewer impacts in each loading cycle. The theory basically is an elaboration of Masri's work which is itself based on the concept of a coefficient of restitution and the conservation of linear momentum. Previous experiments suggest that answers may not be unique in those instances when the theoretically neglected transients cannot be ignored [l].

GENERAL THEORY
The idealized vibroimpact device is shown in Figure 1. Its primary system consists of a linear spring with stiffness K, a viscous dashpot having damping constant C, and mass M excited by the external harmonic load, FO sin Lb. The secondary system is composed of a rigid mass, m, which can move uniaxially in a slot inside the usually much heavier mass h4. The supposed frictionless motion of m is instigated by collisions with M which occur intermittently because of clearance d. Any one impact, typically i, is assumed, reasonably for metals [l, 22,231, to be instantaneous and to be described by a coefficient of restitution, Ri. This coefficient may vary periodically in the general formulation but constant values, RI and R2 corresponding to collisions on the right-and left-hand sides of the slot in Figure 1, will be stipulated in illustrative examples.
Equations of motion and their solutions will be developed in a manner similar to that in reference [l]. The main difference is that the conservation of linear momentum now has to be employed along with the definition of Ri to accommodate the discontinuous velocities at a collision. Displacements, conversely, are invariably continuous and the sinusoidal load sustains motions of the primary and secondary systems which bear similar periodicity conditions to before. However, the periodicity is related more conveniently to the temporal behaviour of the load rather than the displacement of the primary mass. This relationship is illustrated in Figure 2. A motion's periodicity is shown in Figure 2(b) to be To during which time the load in the particular example of Figure 2(a) has undergone two cycles. Meanwhile N impacts between m and M, represented by the spikes in Figure  2(c), have happened. Consequently this motion will be labelled N impacts/two cycles (of the load) with the generalization being obvious. The motion is periodic so that spike YcNtlj at instant t(N+l), equal to (tl + To), iS identical to Y1, YcN+z) at t(N+Z), or (h + To), is identical to Y2, etc. The traditional approach will be taken of considering the periodicity, To, the number of impacts, N, and their occurrence at instants tl, t2, . . . , tN as known. The load's amplitude and phase with respect to the first impact on the other hand will be treated together with the absolute displacement of M at each of the N collisions as unknown. Relationships between these (N + 2) unknowns will be developed next.
where X(t) represents the absolute displacement of the primary mass M. A dot superscript indicates differentiation with respect to time. The solution of equation (1) can be obtained by superimposing free and forced motions [5,241. A typical solution between two consecutive collisions, say the i and the (i + 1) impacts happening at ti and tci+l) respectively, is The temporal origin is taken at tl so that t in equation (2) is the convenient abbreviation of (t -tl). Quantities just after and just before an impact are represented by subscripts a and b respectively. Consequently tia, for example, is the instant immediately after the ith impact. Variables 5, o, r, 7, T and A are given conventionally by Remaining variables are defined as bi =Xi -A sin [ai + T], i=l,2 ,..., N, (a~ is zero because the temporal origin has been taken as the instant of the first impact.) The displacement of M is continuous at collisions so that equation (2) remains valid then. Consequently, the unknown displacement of M at the N collisions within To, Xi, may be obtained by merely substituting their instants of occurrence, ti where i = 1,2, . . . , N, into equation (2). For example, the displacement at the (i + 1)th impact is i=l,2 ,..., N, where, by employing equation (7), The Xci+l) and, from equation (4), A and T are certainly not known in these last three expressions. The (Y~, i = 1, 2, . . . , N, may be calculated from the given driving frequency, 0, and impact instants, ti, by using equation (7). Then the Cl; and Czi may be evaluated by substituting these ai and the determinable parameters of equation (3) into equations (9) and (10). Similar substitutions into equations (5) and (6) leave the final two variables in equation (8), ai and hip requiring Xi, other than the Xi, A and T. The discontinuity in the velocity of A4 at the ith impact may be obtained from the conservation of linear momentum and the definition of Ri.
Linear momentum conservation gives [5,13,231 MXib + &ib = MXia + t&i,, i=1,2 ,..., N, for the ith impact whereas the definition of the coefficient of restitution is The Pib and ii, are the absolute velocities of the secondary mass just before and after the ith impact, respectively. Equations (11) and (12) may be solved straightforwardly to give p=m/M (14) being the (known) mass ratio of the secondary and primary systems. The velocity of m must remain constant between two consecutive impacts because no forces are assumed to act upon it then. Consequently, its velocity ii, between the i and (i + 1) impacts, for example, is simply the ratio of the absolute distance travelled by m and the time elapsed between these two impacts. Mathematically, Zia = (z,i+l,-zi)l(tci+l,-ti), i=1,2 ,..., N. (15) Substituting equation (7) into equation (15) produces

ii, =n[X~i+I,+ Y(i+l)-Xi -Yi]/{a(i+l)-ai),
i=l,2,...,N (16,17) being the relative displacement between m and AI. The Yi is introduced because impacts occur only when the relative displacement is +d/2 for a collision on the right side of M in Figure 1 Consequently, the only unknowns in expression (16)  This last relationship may be substituted into equation (5) to produce an expression for ai in terms of the unknown Xi, A and T. Consequently,  (22) can be substituted into equation (8) to give N linear, simultaneous equations. A perusal of relation (21) should indicate that terms like a(i-1) and a(i+i) have not yet been defined at the two extreme values of i. They have to be obtained from the requirement for periodic motion. It has been seen already that periodicity implies f(N+i) = tl + rr,, t(Nc2) = t2 + TO, . . . , or, in general, if the load undergoes k cycles in time To. Consequently, when i equals N, CY((+~) is simply ((~1 +L?To), or l2T0 in view of equation (7). When i is unity on the other hand, terms like ((Wi -(Yci-i)) in relation (21) may be Written as (a(N+i) -"(N+i-1)) by eII'@Oying equation (25). These last expressions can be evaluated by using the known driving frequency, 0, the assumed periodicity, To, and presumed contact instants ti, i = 1,2, . . . , N. Terms such as Y(N+i-1) can be computed by further presuming the side of it4 at which any particular impact happens. Then equation (18) can be used straightforwardly. All other variables except for the Xi, i = 1, 2, . . . , N, A and T are understood to be given or calculable from equations (3), (14) and (24). Consequently (N + 2) unknowns exist. If trigonometric terms like sin [cz(i+i) +r] in equations (8) and (22) The often cumbersome but calculable coefficients Wriy r = 1, 2, . . . , 6 and i = 1, 2, , . . , N, are given more conveniently in the Appendix. The Xi must still satisfy periodicity condition (24) so that only X1, X2, . . . , XN are actually involved in equation (26). An additional N linear simultaneous equations may be developed by considering the velocity rather than primarily the displacement of M. The velocity of M just before the 5 AND K. McLACHLAN (i +l)th impact, for example, is simply the derivative of equation (2) at t(i+i)b. An alternative expression may be found like that for the velocity just after an impact from linear momentum and the definition of the coefficient of restitution. Relating the two expressions produces the desired equations. From equation (2) with equation (7) and Xci+ijb rather than X(fci+l)b)) being employed. Solving the momentum and restitution equations (11) and (12) for X;.b rather than Xi, gives

J;riib={(Ri-CL)/(l+Ri)}iib+{(l+E.L)/(l+Ri)}ii,,
where p is defined by equation (14). Expressions (16) and (19) for ii, and .& may be substituted into equation (28) to yield . The two resulting equations are independent except in the special, essentially symmetric case when coefficients of restitution are identical and the durations between any three consecutive impacts are the same. This situation is analogous to that observed previously for exclusively rigid impacting mechanisms [l]. Consequently the same procedure as before was adopted: i.e., assume some value for A then calculate the ensuing T. More generally however, the IMSL subroutine LEQTlF was used to solve the four independent equations (26) and (30) numerically by utilizing Gaussian elimination [26]. This course was felt more expedient than solving the equations analytically to produce extremely lengthy expressions. A numerical techt Alternatively, the A and r may have been assumed rather than the (Y~ and the 2N unknowns of Xi and air i = 1,2,. . . , N, calculated, maybe somewhat more precisely, from the 2N equations. However, the resulting equations are linear no longer so that numerical computations are very likely to be much more time consuming.  (26) and (30), g2, was reduced so that E was less than 3% of these values. (The 3% was decreased to 1.5% with virtually no effect in several arbitrarily chosen examples.) All computations were performed throughout the present work with double precision arithmetic on an AMDAHL V7 digital computer. Periodic solutions are grouped conventionally into asymptotically stable or unstable subsets. A periodic motion is asymptotically stable if, after a small perturbation, the original motion is restored eventually [7,13,231. Consequently asymptotic stability would enable machinery performance to be predictable [22,27-291. On the other hand, asymptotic instability implies that a small perturbation will change some aspect of the periodic motion for all time. The asymptotic stability of periodic motions was determined by the Kobrinskii's [23] and Masri [7, 131 for contact surfaces having a single constant coefficient of restitution. Small perturbations at or just after the first impact were related to corresponding changes at or just after the second, third, etc., impacts. This technique was adapted in a companion paper to a rigid impact-pair with different, maybe periodically varying, coefficients of restitution. A similar, somewhat minor extension is possible here. The absolute displacement and velocity of M are no longer prescribed here so they need perturbing along with the previously changeable 2 and 2. Some of the presumed collision instants, and hence the ai, may be modified in accordance with the equations of motion and periodicity conditions. However, these adjustments must be transitional for asymptotically stable motion because the timing sequence must revert eventually to the claimed periodic form. Masri [13] adopted a slightly different approach but with the same end results. He replaced 2, which equation (16) indicates depends upon (yi, by the (Y, of equation (7). Then the initial perturbation and subsequent changes, each denoted by A, were shown to be related by . If all the moduli of the eigenvalues of [P] are strictly less than unity, then a periodic motion is asymptotically stable [13]. Eigenvalues were computed by employing Hessenberg's method in the readily available EIGRF of the IMSL library [26].

WITH PREVIOUS THEORETICAL RESULTS
Computations were performed by using the theory of the last section in order to compare the results with more limited theoretical data available for several particular   A numerical comparison of previous [13] Figure  3 and Table 1 or, regardless of stability, to two unequispaced impacts/cycle in Figure  4. More comprehensive calculations based upon the general theory are presented conventionally in Figure 3 as continuous curves relating values of dK/Fo to the frequency ratio r or 0/w. Particulars of the parameters employed are given at the top of each figure or in Table 1. In addition the periodicity, location, number and distribution of impacts were presupposed from values assumed by previous authors. Consequently, two values of r were considered for the equispaced impact cases of Table 1 [13]. Stability zones displayed in Figure 3 are similar to ones measured for the analogous impact-pair. Associated boundaries are indicated by two curves having a common symbol. Any point within compatible boundaries corresponds to a stable motion with the same periodic form. Further comparisons are possible for two equispaced impacts/cycle alone but they do not add to conclusions from present communications [25]. Table 1 suggests that agreement is normally closest with Masri's results for two unequispaced impacts/cycle. Then the otherwise strong correlations intimate that the single major change in X,,, from 13.79 to 14.77 arises most likely from a printing error in reference [13]. Similar conformity is demonstrated initially for the equispaced impacts but significant differences seem to develop in the eigenvalues. These differences have t Equispaced indicates that the durations between any three consecutive impacts are identical. More general, unequispaced impacts where consecutive durations are unequal will be implied when the distribution of impacts is not stated. The future abbreviation of cycle will infer cyck of the sinusoidal load. 9 been shown [25] to stem largely from a slight slip in former manipulations of the simultaneous equations which affects only the equispaced case. The slip also, plausibly, .explains the premature termination of Masri's equispaced results at r equal to 1.04 on the line d = 37(F0/K) in Figure 3. It does not account on the other hand for Masri's not obtaining the very narrow two unequispaced impacts/cycle stability zone near r equal to 1.00. The maximum eigenvalue modulus is near unity in this zone. Consequently the "all-or-nothing" stability criterion could easily produce opposite conclusions from even minor variations in numerical data. Such variations could well arise from Masri's, apparently, using single rather than, presumably more accurate, double precision arithmetic in numerical computations [7]. Apart from these discrepancies, agreement is reasonable overall in Figure 3 and also at point "a" corresponding to the lone three impacts/cycle example. The absence in this figure of stable equispaced impacts between 0.960 < r < 0.985 supports the previous observation that unevenly spaced impacts are most likely when r is near one. An additional comparison is given in Figure 4 of the results from the degeneration of the general theory into two unequispaced impacts/cycle and those from a Fourier series approach proposed by Sadek [6]. Sadek disregarded the question of stability so that the information is presented without such a determination. Figure 4(a) indicates that values of dK/Fo correlate quite closely whereas large differences occur generally in Figure 4(b) for sin r at frequency ratios, r, near unity. A detailed examination of Sadek's theory indicates that the discontinuity in the velocity of the primary mass at an impact is treated improperly. Consequently, equation (13) of reference [6] implies that the phase angle T is independent of the amplitude F. of the external force. This assertion is refuted clearly by Sadek's own results presented as the solid curves in Figure 4. An independent increase in F. reduces dK/F,, but r remains unaltered at a value, for example, of 0.999. Consequently point B is transposed leftward along the solid curve r = O-999 to position B', say, if the periodic motion, despite a change in the ratio of durations, is still to have two unequispaced impacts/cycle. Points B and B' respectively correspond to points A and A' in Figure 4(b) which obviously relate to different T.
In summary, the theory's conformity with previously meagre results is promising but not convincing. More complete verification is needed. Therefore a mechanical system was built to simulate the mathematical model as far as possible in order to provide a further independent test.

EXPERIMENTAL DETAILS
The experimental analog of the ideal primary system in Figure 1 was composed of items 4, 6 and 18 in Figure 5. Slotted mass 6 with adapters, 4, was supported on two sides by flexible, spring steel strips interposed by a very stiff but light hollow beam, 18. These supports were made so long, 12.00 * 0.08 in, that their ends closest to the slotted mass never moved more than 0*38"* 0.05" from the vertical. Consequently gravitational effects should be negligible. The supports were also much stiffer in torsion than flexure so that the experimental movement was almost uniaxial to conform to conditions imposed by theory. Tightened screws and bolts were used to connect structural components and enable the primary system to move integrally. It is practically impossible to prevent all relative motions, however, so that friction developed at joints. Dissipation was found from conventional free decay and sinusoidal resonance tests [24, 301 to be equivalent to a viscous damping ratio, 6, of 0.0114 f 0.0005. In addition the first and second natural frequencies of the primary system were found to be 19.87 f 0.03 Hz and 430* 2 Hz. Further measurements of select amplitudes and phases of the displacement of the primary mass also indicated good linearity over extreme load characteristics. Hence, the primary system fairly represented a linear oscillator for the frequencies between 10 and 30 Hz employed subsequently. The total weight was equivalent to the tip weights of brackets, slotted mass and accelerometer, labelled 7 in Figure 5, plus 33/140 times the weight of each support for an essentially fundamental cantilevered action [24]. This total was calculated from straightforward measurements to be l-61 lb(f). The corresponding stiffness was approximated for light damping as the square of the measured fundamental radian frequency multiplied by the equivalent weight divided by the gravitational constant. This procedure gave K in Figure 1 as 65.01 lb(f)/in. The secondary system consisted of a dumbbell-like mass, item 5 in Figure 5, strung by two light threads, 10, from a solid frame. It was identical to the mature stainless steel secondary mass described previously in similar impact-pair experiments [l]. The dumbbell reduced extraneous lifting which happened increasingly at higher loading amplitudes but the threads alleviated rotations about the vertical. Hence, the theoretically assumed free, uniaxial motion was duplicated reasonably for the secondary mass which weighed 0,044 f 0.002 lb(f). This mass was effectively rigid too because its lowest natural frequency was 500 times, or more, greater than any loading frequency. High speed photography measurements of coefficients of restitution have been reported earlier [l, 251. Experiments were performed only with magnetically insensitive, mature steel which gave the most consistent results for impact-pairs. As before, unequal coefficients of restitution were achieved by taping one side of the primary mass.
The primary system was connected as in Figure 5 through an impedance head and spring to the electromagnetic shaker. The spring was much weaker than the primary system's effective stiffness so that a force rather than a displacement-like input was obtained. The sinusoidally time varying input, generated to within 2% of the nominal amplitude, was monitored with the aid of the directly coupled impedance head. Compensation for the extraneous mass of this head was achieved by employing the electrical compensation circuit recommended by its manufacturer [31]. Measurements were taken 11 at night to avoid extraneous building vibrations. Absolute displacements and accelerations of the slotted primary mass were observed separately by using a capacitance displacement transducer, 8, and accelerometer, 7. Displacements were determined relative to one of the aluminium brackets, 4. The other bracket maintained symmetry and prevented unbalance. Experimental procedures have been detailed previously [l, 251 and, hence, will not be described again. The only additional limitation was imposed by the restricted force of the shaker [32] which curtailed measurements at the highest periodic impact numbers. Stability had to be demonstrated conclusively again within an arbitrarily imposed but previously satisfactory time limit of 2 min [l].

COMPARISON OF EXPERIMENTAL AND THEORETICAL RESULTS
Experimental and theoretical stability zones are compared in Figures 6-8 with the inverse of dK/Fo in Figure 3 taken as the ordinate. Vertical lines correspond invariably to experimental data whilst the continuous curves present related theoretical results. Coefficients of restitution both average 0.75 in Figures 6 and 7 whereas Figure 8 shows the results of taping one contact surface to produce dissimilar coefficients of 0.75 and 0.61. In actuality, a coefficient of restitution was not constant but decreased generally by about 9% as the frequency ratio, I, increased within the experimental range indicated in the abscissa. Variations from the constant coefficients assumed can be anticipated from analogous impact-pairs to give various but on the whole small fluctuations for the different stability zones considered [ 1, 251. Other parameters observed experimentally and employed in the numerical computations are reported in the figure captions. Figures  6 and 7 relate to the same vibroimpact system; the two impacts/cycle stability zones were separated only to improve clarity. The experimental five impacts/cycle stability zone, distinguished in Figure 7 by vertical lines terminated with dots, was abbreviated by the restricted force capability of the shaker. Corresponding theoretical extremities are indicated simply by dots. Periodicity, impact rate and temporal distribution needed at the start of calculations were obtained from experimental histories like those displayed in Figure 9. (The parts (a)-(c) of Figure 9 correspond to the identically lettered points in Figure 8.) A collision is identified by a sudden acceleration in Figure 9 followed by a gradually decaying envelope as the resulting elastic stress waves diminish. These later transients appeared unimportant providing they were attenuated completely before the subsequent collision [33]. Then the idealization of effectively an acceleration pulse with infinitesimal duration is reasonable. Furthermore, the number and timing of impacts in, for example, Figures 9(a) and (bj can be determined straightforwardly. Difficulties arose, however, when collisions occurred repeatedly and very rapidly on the taped side of the primary mass. This phenomenon is illustrated in Figures 9(c) and (d). It is analogous to the sliding observed for the impact-pair with unequal coefficients of restitution [l]. Such behaviour is outside the scope of the present theory. Figures 6-8 indicate that the non-sliding experimental and theoretical stability zones generally agree well. The largest differences occur in the highest numbered zones and, somewhat surprisingly, in the two equispaced impacts/cycle zone where coefficients of restitution are equal. A close inspection of the equispaced impact zone revealed that the maximum absolute eigenvalue of [P] in equation (32) varied between 0.996 and 0.999. Computations were repeated with single precision arithmetic for this case alone and the ensuing upper stability boundary can be seen from Figure 6(a) to be clearly different. Correlation with experiment is usually poorer with single precision although double precision arithmetic gives an apparently superfluous spike in the practically significant region of r equal to unity. It was impossible, however, to control and load sufficiently finely there to obtain reliable experimental results which might confirm the spike's existence. A similar difficulty was experienced previously [4] which suggests that this point is worthy of further investigation. These experiences together with the narrowness and closeness of all stability zones in the peaks or troughs in Figures 3 and 6-8 indicate that the system's behaviour in the vicinity of r equal to one is very sensitive to small fluctuations in the parameters. Consequently, the traditional two equispaced impacts/cycle representation would appear to have limited applicability in such circumstances. Remaining differences between experiment and theory are similar to those noted previously for an impact-pair [ 11. Consequently, sources of experimental error introduced by the new spring coupling to the shaker and interactions between closer modes of the primary system seem acceptable. Major discrepancies at the largest loads which excite the highest impact rates are largely created again by the uncontrolled lifting of the experimental secondary mass.
Information presented virtually conventionally in Figures 6-8 cannot be applied easily to what is, probably, the most important practical application where vibroimpacts act as a "damper". Then, the secondary mass should oppose and, hence, reduce the originally resonant displacement of the primary system. Attenuation happens only when _&,,,/A, the ratio of the maximum deflections of the primary system with and without impacts, is less than unity. Consequently the alternative representation of Figure 10 FO/(dK) is plotted for a constant r near unity against &,,,/A rather than a variable r, is most appropriate. The figure suggests that the primary system's maximum periodic deflection is attenuated most when r equals unity. It then seems beneficial to have equal coefficients of restitution and two stable equispaced impacts/cycle. The attenuation at a given FO/(dK) diminishes for a small increase in r from unity but, very importantly, amplifications always occur and may be fairly large for a slight decrease in r. Therefore the performance of a traditionally designed, virtually viscous-less impact damper [7, 211 will deteriorate noticeably if r falls consistently below one in practice. Consequently the driving frequency or speed will require close controls with the addition of an impact damper.

CONCLUSIONS
A more general theory has been developed to accommodate any number of repetitive impacts influencing a, formerly, resonant system. Although most incongruences with earlier, sparse predictions of stability may be rationalized, a completely convincing case is not presented for the theory's accuracy. Credibility is enhanced greatly therefore by the close correlation demonstrated between the results from the general theory and comprehensive experiments. The recommended procedure for designing impact dampers is seen to be reasonable for one example of a lightly damped primary system with external sinusoidal loading. This particular example, however, suggests that the operation of the primary system should be maintained at or slightly above its individual fundamental resonance. Careful control of the driving speed may be needed in practice because of the extreme sensitivity displayed by the response in this frequency region.
(A5) The following variables are defined to reduce the length of subsequent expressions:

59
Variables can be computed in the same manner as the coefficients of equations (A2) and (A3) with the additional terms *ia, Zia, A and T determined by the particular periodic solution under consideration. Subscript i in equations (A6) and (A7) varies from 1 to iV.