IImpact Mechanics of Elastic Structures with Point Contact
Róbert Szalai
Department of Engineering Mathematics,University of Bristol, Queen’s Bldg., University Walk,Bristol, BS8 1TR, UK, Email: [email protected] (Dated: January 18, 2014)
Abstract
This paper introduces a modeling framework that is suitable to resolve singularities of impactphenomena encountered in applications. The method involves an exact transformation that turnsthe continuum, often partial differential equation description of the contact problem into a delaydifferential equation. The new form of the physical model highlights the source of singularities andsuggests a simple criterion for regularity. To contrast singular and regular behavior the impactingEuler-Bernoulli and Timoshenko beam models are compared. a r X i v : . [ m a t h . D S ] J a n . INTRODUCTION Impact mechanics is a great concern for engineers, thus many models were developed tounderstand this phenomena [16]. Even simple models of impact [7] lead to complicated pre-dictions, like infinite chatter [11], period adding bifurcations, chaos [1] and non-deterministicmotion [10]. The more elaborate models, however can suffer from convergence problems aseither the time step of the solution decreases [23] or the resolved degrees of freedom isincreased [9]. This signals the need for a better modeling framework for impact phenomena.Most state-of-the-art impact models are finite dimensional and predict infinite contactforces. This happens because at impact the contact point abruptly changes its velocity andhas an infinite acceleration. If the mass of the contact point is not zero, this infinite accel-eration requires an infinite contact force. One example is the standard modal description[5] of linear structures. Each vibration mode is associated with a non-zero modal mass.Truncation of the modal expansion describes the motion as if a finite set of rigid bodies werecoupled with springs and dampers. The arising infinite forces are difficult to handle and arethe source of singularities. In contrast, contact points of elastic bodies have infinitesimallysmall mass. This means that the contact force can stay finite despite infinite accelerations.In this paper we show that finite contact forces are possible under general conditions.In some models, impact is treated as a discrete-time event. Since impacting bodiesmust not overlap, the velocity state of the bodies must be altered at contact, so that theirsubsequent motion avoids overlap momentarily. One such change of velocities can describedby the coefficient of restitution (CoR) model that stipulates that the incident velocity anda rebound velocity of the contact points are opposite and linearly related. This still leaves agreat deal freedom in choosing the rest of the rebound velocities of the structure, thereforethe CoR model must be coupled with additional rules, which can be based on momentumbalances [13, 21] or collocation [22]. The main weakness of the CoR method is that itleads to high-frequency chatter that is proportional to the highest natural frequency of thesystem. This phenomenon restricts the highest natural frequencies that can be included inthe model, because numerical simulation would becomes prohibitively slow. There are a fewextensions to the CoR method that avoid high-frequency chatter [11, 14], however they donot fix the cause of the problem which is the presence of infinite or undefined contact forces.To illustrate this point we compare our method to a CoR model [21] and show how chatter2s eliminated when the contact force is well-defined.Impact can also be modeled by calculating the contact force using impulse response func-tions of the elastic structures at the contact point. This approach is expected to be moreaccurate, however convergence can be as troublesome as for CoR models. For example,Wang and Kim [23] found that in the limit of zero time-step of their algorithm, the contactforce diverges when an Euler-Bernoulli beam impacts an elastic St Venant rod. This diver-gence is related to the approximation of the impulse-response function with finitely manyvibration modes. Similar approaches were used by other authors [4, 6, 24] using time step-ping algorithms and finitely many vibration modes; in some cases the shape of the contactforce was assumed [8], which avoids divergence.In this paper we use a similar technique to impulse-response functions, and transformthe governing equation of the impacting mechanical system to a delay differential equation.The memory term of our equation can be thought of as the convolution integral with theimpulse response function. Depending on the properties of the memory term the model iseither singular or regular. In case of a regular model the finite contact force is uniquelycalculated from a delay-differential equation.
II. MECHANICAL MODEL
We analyze impact mechanics through a converging series expansion of the continuumproblem. We use infinitely many vibration modes x k ( t ) to recover the entire motion of thestructure. Through this expansion the displacement of any point χ of the elastic body canbe written as an infinite sum u ( t, χ ) = ∞ (cid:88) k =1 ψ k ( χ ) x k ( t ) , where ψ k ( χ ) are the mode shapes of the structure [5]. If the motion of the structure isdecoupled into non-resonant modes of vibration, the equation of motion can be written as ¨ x ( t ) + 2 D Ω ˙ x ( t ) + Ω x ( t ) = f e ( t ) + n f c ( t ) , (1)where x = ( x , x , . . . ) T , the mass matrix is assumed to be the identity, Ω = diag( ω , ω , . . . ) and D = diag( D , D , . . . ) , f e ( t ) represents the external force, and f c ( t ) is the contact force.We assume that the natural frequencies scale according to ω k = ω k α , for k (cid:29) . Vector n
3n equation (1) represents the contribution of the modes to the motion of the contact point y ( t ) = n · x ( t ) with n = ( ψ ( χ (cid:63) ) , ψ ( χ (cid:63) ) , . . . ) T , (2)where χ (cid:63) represents the contact point. The method described in this paper is not restrictedto modal equations (1), a more general description can be found in [17]. A. Approximating the contact force
To better understand the impact process we first approximate the contact force assumingthat the impact is infinitesimally short. We assume a single structure that interacts with arigid stop. Contact occurs at t if y ( t ) = 0 . As a first step we calculate a constant contactforce that keeps the stop penetrating the structure after time δt , that is, y ( t + δt ) = 0 , whichforms a boundary value problem. After solving equation (1) the contact force becomes f c = − Cδt α − (cid:0) n · ˙ x ( t − ) (cid:1) , where ˙ x ( t − ) is the vector of modal velocities just before the impact and < C < ∞ is aconstant. When δt → the velocity of the contact point reverses and that corresponds to aunit CoR. The details of the calculation can be found in the Appendix.When evaluating the contact force there are three cases as δt → at the onset of contact.If α < the contact force becomes zero, if α = 1 , the contact force is a finite constant andfor α > the contact force tends to infinity. This simple result implies that for a finitecontact force at least a subsequence of the natural frequencies ω k must scale at most linearlyas k → ∞ .For a system composed of an elastic body which strikes a rigid stop, the impact shouldchange the momentum of the elastic body. Our calculation shows that the change of mo-mentum of the elastic body is zero, i.e., lim δt → ( f c δt ) = 0 as the contact time tends to zero.This implies that the impact must occur during a non-zero and finitely long time-interval. III. MODEL TRANSFORMATION
To accurately calculate the contact force as a function of time in the continuum problemwe transform the infinite dimensional system (1) into a delay equation. Delay terms can nat-urally arise from traveling wave solutions of partial differential equations [15, 18]. However4ispersion might prevent one to write down such solutions. Instead we use the Mori-Zwanzigformalism as is described for mechanical systems in [17] and obtain a time-delay model.Our aim is to find a self-contained equation that exactly describes the evolution of y ( t ) =( n · x ( t ) , n · ˙ x ( t )) T . We call y the vector of resolved variables. The first step in the processis to transform (1) into a first-order form ˙ z ( t ) = Rz ( t ) + n f c ( t ) + f ( t ) , R = I − Ω − D Ω , (3)where f ( t ) = ( , f e ( t )) T . Note that R can also represent any convergent expansion ofthe continuum problem, including numerical schemes such as finite difference methods and y can include any finite number of variables [17]. To arrive at a model that describes theevolution of the resolved variables y , we construct a projection with a finite dimensionalrange with the help of the matrices V = n T n T and W = m m , where the vector m is chosen such that m · n = 1 and its components obey [ m ] j = 0 for j > M < ∞ . To simplify our analysis we also assume that the columns of W span aninvariant subspace of R . The resolved variables now can be expressed as y = V z , andthe projection and the reciprocal projection matrices become S = W V and Q = I − S ,respectively. Further, we assume that the initial condition of (3) is specified at t = 0 andthat there is no contact force initially. According to [17], with this notation the governingequation for the resolved variables becomes dd t y ( t ) = Ay ( t ) + L ∞ f c ( t ) + ˆ t d τ L ( τ ) dd t [ f c ( t − τ )] + g ( t ) , (4)where A = V RW , the memory kernel is a function of bounded variation, L ( τ ) = ˆ τ (cid:0) V e RQ θ (0 , n ) T − L ∞ (cid:1) d θ, L ∞ = AV R − (0 , n ) T , and the forcing term is g ( t ) = V RQ e R t ( x (0) , ˙ x (0)) T + ˆ t ( V R − AV ) e R τ f ( t − τ )d τ. The integral in (4) is meant in the Riemann-Stieltjes sense. This means that discontinuitiesof L ( τ ) at τ i represent discrete values of dd t f c ( t − τ i ) .5 V. CONTACT FORCE CALCULATION
In what follows, we interpret the meaning of (4) for impact problems. In the simple casewhen impact occurs with a stationary stop the contact point must have both a constantposition y = y and zero velocity y = 0 . This means that at the time of initial contactthe acceleration of the contact point becomes infinite as the velocity resets to zero. Despitethe infinite acceleration, the zero mass of the contact point guarantees a finite contact force.We show that in the transformed model the zero mass of contact point is equivalent to thecondition lim τ → [ L ( τ )] (cid:54) = 0 . (5)In other words, if condition (5) is satisfied the contact force is finite. This is a similar, butmore general condition that has been derived by Wang and Kim [23]. For multiple resolvedcoordinates equation (5) must hold for the coordinates that experience a discontinuity atimpact, e.g., the velocities. To show that our conclusion holds we integrate equation (4)through the initial contact to get y ( t +0 ) − y ( t − ) = L + (cid:0) f c ( t +0 ) − f c ( t − ) (cid:1) , where L + = lim τ → L ( τ ) and t − signals a limit from the left and t +0 a limit from the rightof the impact. Because the contact force before the impact f c ( t − ) = 0 and the velocity ofthe contact point after the impact y ( t +0 ) = 0 , it follows that the initial contact force is f c ( t +0 ) = − y ( t − ) (cid:2) L + (cid:3) . (6)During contact y ( t ) = y = ( y , T is constant, which can be substituted into (4) to findout the contact force. Rearranging the resulting equation yields (cid:2) L + (cid:3) dd t f c ( t ) = − [ L ∞ f c ( t ) + A ¯ y + g ( t )] − ˆ t d τ [ L ( τ )] dd t [ f c ( t − τ )] . (7)Equation (7) is a delay-differential equation that contains the history of f c ( t ) , which meansthat previous impacts have a great influence on the evolution of the contact force. Clearly,the contact force is not a continuos function, therefore the derivative of its history canbecome infinite. A short calculation shows that a jump of magnitude f jump c of f c at t − τ (cid:63) contributes a finite value d / d τ L ( τ (cid:63) ) f jump c to the integral in (7).6n conservative systems where shock waves are present, e.g., for an undamped string, L ( τ ) has further isolated discontinuities. Let τ d (cid:54) = 0 be the position of such a discontinuityof L ( τ ) . If at time t = t + τ d the two bodies are in contact the contact force develops afurther jump. This can be seen by integrating equation (7) for the infinitesimal time interval [ t − , t +1 ] . The integral of the integral on the right side of (7) becomes ˆ t +1 t − ˆ t d τ [ L ( τ )] dd t [ f c ( t − τ )] d t = (cid:20) ˆ t d τ [ L ( τ )] f c ( t − τ ) (cid:21) t +1 t − , (8)while the other terms are continuous and their integral vanishes. The right side of equation(8) is regular because all of its terms are finite. The discontinuity of L ( τ ) at τ d contributes (cid:0) L ( τ + d ) − L ( τ − d ) (cid:1) f c ( t − τ d ) to the integral (8). Note that t − τ d = t , therefore (8) evaluatesto (cid:0) L ( τ + d ) − L ( τ − d ) (cid:1) f c ( t +0 ) (as f c ( t − ) = 0 ), which means that the contact force at t alsodevelops a further discontinuity L + (cid:0) f c ( t +1 ) − f c ( t − ) (cid:1) = (cid:0) L ( τ + d ) − L ( τ − d ) (cid:1) f c ( t +0 ) . A. Numerical solution of the reduced model
The non-smooth delay-differential equations (4,6,7) are somewhat unusual and thereforestandard numerical techniques are not directly applicable to them. In what follows, we use asimple explicit Euler method and the rectangle rule of numerical integration to approximatethe solution of (4) for cases when (cid:2) L + (cid:3) (cid:54) = 0 . We assume that time is quantized in ε chunks, so that y q = y ( qε ) , f c,q = f c ( qε ) , where q = 0 , , , . . . . If there is no contactand consequently no contact force ( f c,q = 0 ), the only unknown is the state variable y q .Therefore the evolution of the resolved variables is given by y q +1 = y q + ε (cid:0) Ay q + g ( qε ) (cid:1) + q − (cid:88) j =0 ( L j +1 − L j ) ( f c,q − j − f c,q − j − ) , (9)where L j = L ( jε ) and L = L + . Contact of the impacting bodies is detected, when (cid:2) y q +1 (cid:3) ≤ y . In this case the resolved variables are kept constant with y q +1 = y . Also,equation (6) is applied at the onset of contact, so that the initial contact force becomes f c,q +1 = − (cid:2) y q (cid:3) (cid:2) L + (cid:3) . (10)7 ( t , ) u ( t ,1) Figure 1: Impacting cantilever beam.
The subsequent values of the contact force are calculated by f c,q +1 = f c,q − ε (cid:2) L + (cid:3) [ L ∞ f c,q + A ¯ y + g ( qε )] − (cid:2) L + (cid:3) q − (cid:88) j =0 [ L j +1 − L j ] ( f c,q − j − f c,q − j − ) . (11)which is the discretized counterpart of (7). If f c,q +1 as predicted by equation (11) becomesnegative, we set f c,q +1 = 0 and continue the calculation with (9). V. IMPACTING CANTILEVER BEAM MODELS
Our theory sets a criterion in the form of equation (5) for the regularity of the mechanicalmodel, which can be used to test different models of elastic structures. We consider theexample of a cantilever beam in Fig. 1 described by two different models. Through ourcalculation it becomes clear why the Euler-Bernoulli model often used in impact models[6, 9, 24] exhibits signs of singularity.As the underlying elastic structure, consider the Euler-Bernoulli cantilever beam model ∂ u∂t = − ∂ u∂ξ + ψ f e ( t ) , u ( t,
0) = ∂u ( t, ξ ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =0 = ∂ u ( t, ξ ) ∂ξ (cid:12)(cid:12)(cid:12)(cid:12) ξ =1 = 0 , (12)with ∂ u ( t,ξ ) / ∂ξ | ξ =1 = f c ( t ) , where u ( t, ξ ) represents the deflection of the beam. The naturalfrequencies of (12) are determined by the equation √ ω k cosh √ ω k = 0 , while the modeshape values at the end of the beam are described by n = (2 , − , , − , . . . ) T . (13)On the other hand the Timoshenko beam model is represented by ∂ u∂t = βγ (cid:18) ∂ u∂ξ − ∂φ∂ξ (cid:19) + ψ f e ( t ) ,∂ φ∂t = β ∂ φ∂ξ + β γ (cid:18) ∂u∂ξ − φ (cid:19) ,u ( t,
0) = φ ( t,
0) = φ ( t,
1) = 0 , (14)8 uler-BernoulliTimoshenko Euler-Bernoulli Timoshenko (a) (b) -0.15-0.1-0.0500.150.10.050 0.5 1 1.5 2 2.5 3 Figure 2: (color online) Graph of [ L ( τ )] for the Euler-Bernoulli and the Timoshenko beam models.The inset shows that the Timoshenko beam model converges to a function where (cid:2) L + (cid:3) (cid:54) = 0 , whilethe Euler-Bernoulli model is singular since its L ( τ ) is continuos. The parameters are β = 4800 , γ = / . and ∂ / ∂ξ u ( t, | ξ =1 − φ ( t,
1) = f c ( t ) , where f e ( t ) is an external forcing through the secondmode shape ψ , and φ is the rotation angle of the cross-section of the beam. The resolvedcoordinates in both cases are y = u ( t, and y = ∂ / ∂t u ( t, . In case of the Timoshenkobeam, modal decomposition in the form of (1) is not practical, instead we use Chebyshevcollocation [19] to discretize the system with N number of collocation points. This is possible,since our formulation is not restricted to modal decomposition, the matrix R can representany form of discretization. The governing equations (12,14) are conservative, thus to obtaina decaying solution we add modal damping ratios D k = / to both systems (12,14) afterbeing discretized. The result of our calculation is shown in Fig. 2. It can be seen that forthe Euler-Bernoulli model [ L ( τ )] is continuos while for the Timoshenko model [ L ( τ )] isdiscontinuous. This result is in accordance with the fact that the scaling exponent of thehighest natural frequencies for the Euler-Bernoulli model is α = 2 and for the Timoshenkomodel α = 1 [20]. This means that the Euler-Bernoulli model is not suitable for impactcalculations. 9 a) (c)(d)(b) Figure 3: (color online) Vibrations of the impacting Timoshenko beam using two impact models.The rigid stop as illustrated in Fig. 1 is placed at y = − . , and the forcing is f e ( t ) = 30 cos (13 t ) through the second mode. Trajectories of the reduced model (4,6,7) are shown in dark (blue) andthe solution of the CoR model (1,15) is represented by light (green) lines for comparison. The timestep used to solve the reduced model is ε = 3 . × − and the number of collocation points used tosolve the CoR model is N = 20 . Panels (a,b) show that the solution converges to a periodic orbit.A single period of the solution is illustrated in panels (c,d). A. The Timoshenko model
In the previous section we have shown that the Timoshenko beam model (11) satisfiesour criterion (5), which guarantees that equations (4,6,7) are well defined and that thecontact force is finite during impact. To confirm that this is indeed the case we simulate theimpacting cantilever beam with the Timoshenko beam model (14) as illustrated in figure1. The rigid stop is placed at y = − . , so that the beam contacts the stop whenthe position of its tip reaches u ( t,
1) = − . . In all our simulations we use the initialconditions u (0 , ξ ) = ∂ / ∂t | u ( t, ξ ) | t =0 = 1 . ψ ( ξ ) , where ψ is the first mode shape of thestructure and it is normalized by ψ (1) = 1 . We also force the beam through its second modewith f e ( t ) = 30 cos (13 t ) . The numerical solution is obtained using equations (9,10,11). Inaddition to our method we also simulate the dynamics using the CoR model described in[21]. This comparison highlights that high-frequency chatter is eliminated when our methodis used. 10 a)(b)(c)(d)(e) Figure 4: (color online) A sequence of two impacts within a period of the periodic solution in Fig.3. The solution of the CoR model (1,15) as shown by the light (green) lines is highly oscillatory.The dark (blue) lines show that the reduced model (4,6,7) is similarly accurate and eliminates thehigh-frequency chatter of the CoR model. (a,b). Insets (d,e) show the small-scale dynamics of theCoR model. The contact force is finite and continuos after the initial contact (c).
The CoR model that we are using for comparison describes the impact as an infinitesi-mally short process. During the impact an impulse is applied at the contact point that altersthe velocity state of the body. The magnitude of this impulse is determined by the desiredrebound velocity of the contact point which is − C R times the incident velocity. Assumingthat the equation of motion is in the form of (1) and impact occurs at t , the after-impactvelocity of the structure is ˙ x ( t +0 ) = (cid:18) I − (1 + C R ) n ⊗ nn (cid:19) ˙ x ( t − ) , (15)where ⊗ means the outer product between vectors. The position of the structure remainsthe same throughout the impact: x ( t +0 ) = x ( t − ) . The CoR model is physically questionablebecause it treats the impact as an infinitesimally short event. It can however, reproducemost experimental observations [12]. The source of high-frequency chatter can be explainedby equation (15), which stipulates that the change in modal velocities is proportional to a11onstant times vector n . Elements of n corresponding to high frequency modes are of thesame magnitude as for low frequency modes (e.g., see equation (13)), which means that afteran impact the tip of the beam acquires a high frequency vibration. Due to this vibrationanother impact is likely to follow shortly and repeatedly, which results in chatter that hasroughly the same frequency as the highest vibration mode of the structure. This is illustratedby the light (green) lines in figure 4.In contrast to the CoR method (1,15) our method as solved by equations (9,10,11) elim-inates chatter and produces a much smoother result, which are shown by dark (blue) linesFigs. 3 and 4. On the larger scale in figure 3 the two solutions roughly coincide, while onthe smaller scale in figure 4 the high frequency chatter is apparent and would increase infrequency if more vibration modes or collocation points on the beam were used. This highfrequency chatter can stall numerical simulations, while the time-step in our method is notaffected by the inclusion of higher natural frequencies.The contact force in the CoR model is infinite at times of contact and hence it cannotbe calculated. On the other hand our method allows the calculation of the contact force,which is finite as shown in figure 4(c). The contact force even becomes a smooth functionof time after the onset of contact. VI. CONCLUSIONS
In this paper we introduced a new way of modeling the impact mechanics of elasticstructures. With our method regularity of the model can be predicted and a finite andpiecewise continuos contact force can be calculated. The key to this result is that thedelay equation description preserves the infinite dimensional nature of the mechanics andthe zero mass of the contact point. The results presented in this paper open a significantnumber of new avenues of research. Models that show non-deterministic behavior such as thePainleve paradox [10] might be regularized through our method. The strong dependence ofdynamical phenomena on the number of underlying dimensions [3] could also be eliminated,since our framework considers all the infinite dimensions. Further, the bifurcation theoryof non-smooth delay equations requires attention in order to understand the implicationsof our regularized impact mechanics, especially in how far it is an improvement over finitedimensional models. 12 cknowledgments
The author thanks Gábor Stépán, who brought his attention to the work of Chorin etal. [2]. He also thanks Alan R. Champneys, John Hogan and Thibaut Putelat for usefuldiscussion and comments on the manuscript. Corrections to the text by Galit Szalai aregreatly appreciated.
Appendix A: Contact force asymptotics
In this appendix we approximate the contact force at the onset of an impact. We assumethat an impact takes place at t = t . To help the notation we define (cid:3) − = (cid:3) ( t ) and (cid:3) + = (cid:3) ( t + δt ) , where (cid:3) stands for any dependent variable. We aim to calculate a constantcontact force f c that allows the elastic body to overlap with the rigid stop for an exactly δt long time interval. As δt tends to zero the overlap is removed, hence the calculated f c forcetends to the actual contact force.To calculate this constant f c one needs to solve n · x + for f c . Due to the linearity of equation (1), the motion depends linearly on the contact force.Therefore expanding the constraint n · x + at f c = 0 we get the exact equation n · x + f c =0 + f c N (cid:88) k =1 ψ k ( χ (cid:63) ) ∂x + k ∂f c . (A1)The derivatives in (A1) are calculated in closed form as ∂x + k ∂f c = − ψ k ( χ (cid:63) ) ω k (cid:32) e − D k ω k δt (cid:112) − D k sin (cid:18) ω k (cid:113) − D k δt (cid:19) + e − D k ω k δt cos (cid:18) ω k (cid:113) − D k δt (cid:19) − (cid:33) . (A2)When calculating the derivatives (A2) for δt = 10 − , we get a vanishing sequence as isillustrated in Fig. 5(a).We assume that before the impact, the structure has a smooth motion so that the dis-placement without the contact force can be approximated as n · x + f c =0 ≈ δt n · ˙ x − . a) (b) Figure 5: (a) Coefficients (A2) for δt = 10 − , D k = 0 and ω k = (cid:0) kπ − π (cid:1) . (b) The sum of theinfinite series of coefficients (A2) is shown by the continuos (blue) line. The dashed (red) line is theestimate of N , in Eqn. (A4) that overestimates the sum. To quantify how the derivatives in Eqn. (A2) vanish we asymptotically expanded them as ∂x + k ∂f c = 12 ψ k ( χ (cid:63) ) δt (cid:18) − D k ω k δt − − D k ω k δt + · · · (cid:19) , Substituting these two estimates into (A1) we get the contact force f c ≈ − n · ˙ x − δt (cid:80) Nk =1 ψ k ( χ (cid:63) ) = − n · ˙ x − N δtψ k ( χ (cid:63) ) , (A3)up to the leading order, where the limit N of the summation equals the smallest k for which(A2) vanishes, and ψ k ( χ (cid:63) ) is the average value of ψ k ( χ (cid:63) ) , for k ≤ N . N can be calculatedas the smallest k such that the expanded coefficients become small δt ∂x + k ∂f c ≈ − D k ω k δt − − D k ω k δt < η (cid:28) . Using the asymptotic scaling of the natural frequencies ω k = ω k α and assuming zero damp-ing ( D k = 0 ) we get N > (cid:18) √ − ηω (cid:19) α δt − α , (A4)therefore the restoring force is f c = − Cδt α − n · ˙ x − , where C ≈ (cid:16) √ − ηω (cid:17) − α ψ k ( χ (cid:63) ) − . In order to check validity of our estimate of N we plottedthe sum of derivatives in (A2) and compared to the estimate (A4) in Fig. 5(b).14s the last step we calculate the change in modal velocities ˙ x + k = ˙ x − k − Cδt α − n · ˙ x − ∂ ˙ x + k ∂f c ≈ ˙ x − k − ψ k ( χ (cid:63) ) Cδt α n · ˙ x − , which means that if δt → , there is no change in individual modal velocities, except when α → ∞ , that is the case of a rigid body. However, when calculating the velocity of theimpacting point after the impact we have n · ˙ x + = − n · ˙ x − . [1] W. Chin, E. Ott, H. E. Nusse, and C. Grebogi. Grazing bifurcations in impact oscillators. Phys. Rev. E , 50(6):4427–4444, 1994.[2] A. J. Chorin, O. H. Hald, and R. Kupferman. Optimal prediction and the mori-zwanzigrepresentation of irreversible processes.
Proc. Natl. Acad. Sci. U. S. A. , 97(7):2968–2973,2000.[3] M. di Bernardo, C. Budd, A. R. Champneys, and P. Kowalczyk.
Piecewise-smooth DynamicalSystems: Theory and Applications . Springer, 2008.[4] G. R. Evans, B. C. Jones, A. J. McMillan, and M. I. Darby. A new numerical-method for thecalculation of impact forces.
J. Phys. D-Appl. Phys. , 24(6):854–858, 1991.[5] D. J. Ewins.
Modal Testing: Theory, Practice and Application (Mechanical Engineering Re-search Studies: Engineering Dynamics Series) . Wiley-Blackwell, 2000.[6] A. Fathi and N. Popplewell. Improved approximations for a beam impacting a stop.
J. SoundVibr. , 170(3):365–375, 1994.[7] Y. A Khulief and A. A. Shabana. A continuous force model for the impact analysis of flexiblemultibody systems.
Mech. Mach. Theory , 22(3):213–224, 1987.[8] R. S. Langley. The analysis of impact forces in randomly vibrating elastic systems.
J. SoundVibr. , 331(16):3738–3750, 2012.[9] J. Melcher, A. R. Champneys, and D. J. Wagg. The impacting cantilever: modal non-convergence and the importance of stiffness matching.
Philosophical Transactions of the RoyalSociety A: Mathematical, Physical and Engineering Sciences , 371(1993), 2013.
10] A. Nordmark, H. Dankowicz, and A. R. Champneys. Friction-induced reverse chatter in rigid-body mechanisms with impacts.
IMA J. Appl. Math. , 76:85–119, 2011.[11] A. B. Nordmark and P. T. Piiroinen. Simulation and stability analysis of impacting systemswith complete chattering.
Nonlinear Dyn. , 58(1-2):85–106, 2009.[12] M. Oestreich, N. Hinrichs, and K. Popp. Dynamics of oscillators with impact and friction.
Chaos, Solitons & Fractals , 8(4):535–558, 1997.[13] H. Palas, W. C. Hsu, and A. A. Shabana. On the use of momentum balance and the assumedmodes method in transverse impact problems.
J. Vib. Acoust.-Trans. ASME , 114(3):364–373,1992.[14] D. J. Segalman, A. M. Roy, and M. J. Starr. Modal analysis to accommodate slap in linearstructures.
J. Vib. Acoust.-Trans. ASME , 128(3):303–317, 2006.[15] G. Stépán and Zs. Szabó. Impact induced internal fatigue cracks. In
Proceedings of theDETC’99 , 1999.[16] W. J. Stronge.
Impact Mechanics . Cambridge University Press, 2000.[17] R. Szalai. Modelling elastic structures with strong nonlinearities with application to stick-slipfriction.
Proc. R. Soc. A , 470(2161), 2014.[18] D. Takács, G. Orosz, and G. Stépán. Delay effects in shimmy dynamics of wheels with stretchedstring-like tyres.
Eur. J. Mech. A-Solids , 28(3):516–525, 2009.[19] L. N. Trefethen.
Spectral Methods in MATLAB . SIAM Philadelphia, 2000.[20] N. F. J. van Rensburg and A. J. van der Merwe. Natural frequencies and modes of a timoshenkobeam.
Wave Motion , 44(1):58–69, 2006.[21] C. P. Vyasarayani, J. McPhee, and S. Birkett. Modeling Impacts Between a Continuous Systemand a Rigid Obstacle Using Coefficient of Restitution.
J. Appl. Mech.-Trans. ASME , 77(2),2010.[22] D. J. Wagg and S. R. Bishop. Application of non-smooth modelling techniques to the dynamicsof a flexible impacting beam.
J. Sound Vibr. , 256(5):803–820, 2002.[23] C. Wang and J. Kim. New analysis method for a thin beam impacting against a stop basedon the full continuous model.
J. Sound Vibr. , 191(5):809–823, 1996.[24] X. C. Yin, Y. Qin, and H. Zou. Transient responses of repeated impact of a beam against astop.
Int. J. Solids Struct. , 44(22-23):7323–7339, 2007., 44(22-23):7323–7339, 2007.