# Nonsmooth Modeling of Vibro-Impacting Euler-Bernoulli Beam

Authors: C. P. Vyasarayani; Sukhpreet Singh Sandhu; John McPhee
Journal: Advances in Acoustics and Vibration
 (2012)

---

## Abstract

A new technique to simulate nonsmooth motions occurring in vibro-impacting continuous systems is proposed. Sticking motions that are encountered during vibro-impact simulation are imposed exactly using a Lagrange multiplier, which represents the normal reaction force between the continuous system and the obstacle. The expression for the Lagrange multiplier is developed in closed form. The developed theory is demonstrated by numerically simulating the forced response of a pinned-pinned beam impacting a point-like rigid obstacle.


---

## Body

## 1. Introduction

The dynamics of continuous systems like beams, strings, and rods subjected to impacts caused by motion-limiting constraints have been studied extensively. A representative set of studies can be found in [1–17]. The impact of these systems is modelled using two different methods.(1)
Force integration method (penalty method) [2, 4–6, 10–14]: in this method, a contact force proportional to the penetration of the beam at the contact point is applied at the contact location. This method can be used to model linear or nonlinear stiffness in the obstacles.(2)
Nonsmooth method (coefficient of restitution (CoR) approach): this approach treats the obstacle as rigid and uses the CoR assumption to relate the states before and after impacting neglecting the change in configuration during the contact. A modal form of CoR rule for continuous systems was proposed by Wagg and Bishop [1, 7–9] and further developed by Vyasarayani et al. [18] to make the formulation energy-consistent.The penalty method is best suited for simulating impacts that are compliant but, when used to simulate rigid impacts, they suffer from long simulation times due to the abrupt discontinuity in stiffness arising at the event of impact [19], requiring special numerical treatment [20, 21]. The COR method is an attractive method for simulating the impact with rigid point-like obstacles due to its several advantages. The modeling process is simple because when the impact occurs one simply uses the CoR rule to obtain the post-impact modal states from the pre-impact modal states. An energy-consistent modal form of CoR rule has been proposed by Vyasarayani et al. [18] and was shown to perform better in terms of simulation time than the penalty approach. In the work of Vyasarayani et al. [18], the authors have not discussed the modeling of sticking motions.The aim of this work is to highlight and address the issues associated with modeling sticking motions that can occur when a continuous system impacts a rigid obstacle. During the impact of a continuous system with a rigid obstacle, before a sticking motion occurs, it is preceded by a chatter sequence. In a chatter sequence [22] the beam impacts the obstacle repeatedly, and the time between the two successive impacts (Δτi) continuously decreases. Simulating the chatter sequence without assuming sticking when Δτi≪1 can be computationally prohibitive. For lumped parameter systems, Cusumano and Bai [23] proposed to monitor Δτi in a chatter sequence and, if this time falls below a certain tolerance, the structure is assumed to stick with the obstacle. The reaction at the interface of the stuck mass and obstacle is monitored and when it crosses zero the contact phase is assumed to be over and the mass is released. The sticking motion can also be interpreted as motions occurring with zero coefficient of restitution in the impact model. However, for continuous systems the method proposed by Cusumano and Bai [23] cannot be directly applied.Sticking motions in continuous systems are usually simulated by using mode-switching method [24, 25], wherein when the beam is assumed to be stuck to the obstacle the basis functions of the nonimpacting motion are changed to a different basis, usually the mode shapes of the beam with a pin at the location of the obstacle. The reaction at the pin is monitored to see when it crosses through zero and then the basis is again changed from contact phase to noncontact phase. In this paper, we propose a new method to simulate the beam motion during the sticking phase which does not require the reevaluation of mode shapes with a pin at the location of impact, thus enabling the efficient handling of sticking motions within the CoR method. This is done by deriving closed form expressions for the reaction force at the pin such that the gap constraint is exactly enforced. A suitable numerical example will be shown to demonstrate this method on a vibroimpacting simply supported beam.


## 2. Mathematical Modeling

We consider a simply supported beam with a motion-limiting constraint. The beam-stop system is shown in Figure1. The beam is subjected to a uniformly distributed harmonic force F*sin(ω*t). The stop that limits the motion of the beam is located at Xi. The governing partial differential equation (PDE) for the beam using Euler-Bernoulli beam theory, without considering the impact event, is as follows:
(1)EI∂4Z∂X4+C∂Z∂t+ρA∂2Z∂t2=F*sin⁡(ω*t),Z(Xi,t)≤D,
where E is the Young’s modulus, I is the area moment of inertia, ρ is the mass density, A is the cross-sectional area, C is the viscous damping coefficient per unit length, and t is time. The pinned-pinned boundary conditions are enforced as follows:
(2)Z(0,t)=Z(L,t)=0,∂2Z(0,t)∂X2=∂2Z(L,t)∂X2=0.
We scale the equation of motion by introducing the following parameters:
(3)x=XL,d=DL,w=ZL,τ=tEIρAL4,ω=ω*ρAL4EI,F=F*L3EI,c=CLEIρAL4.
The equation of motion in scaled form is as follows:
(4)∂4w∂x4+c∂w∂t+∂2w∂τ2=Fsin⁡(ωτ),w(xi,τ)≤d.
The boundary conditions are
(5)w(0,τ)=w(1,τ)=0,∂2w(0,τ)∂x2=∂2w(1,τ)∂x2=0.
The initial conditions are as follows:
(6)w(x,0)=w0(x),w˙(x,0)=w˙0(x).
Using the expansion theorem, a solution to (4) is assumed to be of the following form:
(7)w(x,τ)=∑j=1∞‍Wj(x)qj(τ),
where Wj(x)=Csin(jπx) and qj(τ) correspond to the jth undamped mode shape and the corresponding modal coordinate for the beam, respectively. Substituting the above form of solution into (4), we get
(8)∑j=1∞‍d4Wj(x)dx4qj(τ)+∑j=1∞‍Wj(x)q˙j(τ)+∑j=1∞‍Wj(x)q¨j(τ)=Fsin⁡(ωτ).
The above equation is multiplied by Wk(x) and then integrated over the spatial domain, which yields the following expression:
(9)∑j=1∞‍qj(τ)∫01d4Wj(x)dx4Wk(x)dx+∑j=1∞‍q˙j(τ)∫01Wj(x)Wk(x)dx,∑j=1∞‍q¨j(τ)∫01Wj(x)Wk(x)dx=Fsin⁡(ωτ)∫01Wk(x)dx,k=1,2,…,∞.
We now evaluate the following integral:
(10)∫01Wj(x)Wk(x)dx=C2∫01sin⁡(jπx)sin⁡(kπx)dx=C22δjk,
where δjk is Kronecker’s delta, defined as:
(11)δjk={1,ifj=k,0,ifj≠k.
Since the amplitudes of the mode shapes are arbitrary, we choose the amplitude of mode shapes (C) to be 2. This procedure of choosing the amplitudes of mode shapes such that the right-hand side of (10) that becomes δjk is called mass-normalization. By using the mass-normalized mode shapes, it can be shown that:
(12)∫01d4Wj(x)dx4Wk(x)dx=ωj2δjk=(jπ)4δjk.
The coupled ordinary differential equations (ODEs) given by (9) can be decoupled using the orthogonality relations ((10) and (12)). The uncoupled ODEs, after truncating the system to N modes, are as follows:
(13)q¨j(τ)+2ζjωjq˙j+ωj2qj(τ)=αjFsin⁡(ωτ),j=1,2,…,N,
where αj=∫01Wk(x)dx and ζj=c/(2ωj). To obtain the modal initial conditions, we substitute the series solution given by (7) into (6), which gives the following relations for displacement and velocity initial conditions:
(14)w(x,0)=∑j=1∞‍Wj(x)qj(0),w˙(x,0)=∑j=1∞‍Wj(x)q˙j(0).
Multiplying both sides of (14) by Wk(x), and then integrating over the domain, results in the following modal initial conditions:
(15)qj(0)=∫01w(x,0)Wj(x)dx,q˙j(0)=∫01w˙(x,0)Wj(x)dx.Figure 1
Schematic of the physical system.


### 2.1. Impact Modeling

Equation (13) when integrated with appropriate initial conditions can simulate the beam motion until the impact. Let τ- be the time just before the impact (pre-impact time). Pre-impact states, w(x,τ-) and w˙(x,τ-), are known. We have to find the states immediately after impact (post-impact), w(x,τ+) and w˙(x,τ+). Using the classical Newton CoR assumption for rigid impacts, the configuration of the beam is assumed not to change during impact [1] (as rigid impacts are assumed to occur in an infinitesimal time), we have
(16)w(x,τ+)=w(x,τ-).
The post-impact modal displacement initial conditions are
(17)qj(τ+)=qj(τ-),j=1,2,…,N.
We also know that the beam velocity at the location of impact is altered according to the classical CoR rule:
(18)w˙(xi,τ+)=-Rw˙(xi,τ-),
where R is the coefficient of restitution. We must now obtain the post-impact velocity distribution, w˙(x,τ+), satisfying (18). The procedure for obtaining the post-impact velocity distribution is explained in detail in Vyasarayani et al. [18]; we present a brief explanation here. Let us assume that at the instant of impact (τ-) an impulse of magnitude P acts on the beam at the impact location. The equation of motion of the beam can be written as:
(19)∂4w∂x4+c∂w∂t+∂2w∂τ2=Pδd(x-xi)δd(τ-τ-).
We assume the same pinned-pinned boundary conditions for the above PDE. Substituting the series solution given by (7) into the above PDE, we have
(20)∑j=1N‍d4Wj(x)dx4qj(τ)+∑j=1N‍Wj(x)q˙j(τ)+∑j=1N‍Wj(x)q¨j(τ)=Pδd(x-xi)δd(τ-τ-).
Integrating the above equation between τ- and τ+, we get
(21)∑j=1N‍d4Wj(x)dx4∫τ-τ+qj(τ)dτ+∑j=1N‍Wj(x)∫τ-τ+q˙j(τ)dτ+∑j=1N‍Wj(x)∫τ-τ+q¨j(τ)dτ=Pδd(x-xi)∫τ-τ+δd(τ-τ-)dτ.
The post-impact displacement distribution will remain the same as the pre-impact distribution, as the impact is assumed to occur in an infinitesimal amount of time, that is, w(x,τ+)=w(x,τ-). Hence we have qj(τ+)=qj(τ-) and thus the first integral in (21) becomes zero. Integrating the second term with respect to time, and from the definition of Dirac’s delta function, ∫τ-τ+δd(τ-τ-)dτ=1, (21) becomes
(22)∑j=1NWj(x)(q˙j(τ+)-q˙j(τ-))dt=Pδd(x-xi).
Now substituting for w˙(x,τ)=∑j=1N‍Wj(x)q˙j(τ), we get:
(23)w˙(x,τ+)=w˙(x,τ-)+Pδd(x-xi).
From the above equation we can say that the effect of applying an impulse on the Euler-Bernoulli beam results in a sudden change in velocity of the beam, immediately after the application of impact. The resulting velocity distribution after impact is nothing but a Dirac’s function in space added to the pre-impact velocity distribution. The scaling parameter P (magnitude of impulse) in (23) is to be found such that the post-impact velocity at the impact location is as described by (18). Multiplying both sides of (23) by Wk(x) and integrating over the domain results in:
(24)∑j=1N‍q˙j(τ+)∫01Wj(x)Wk(x)dx=∑j=1N‍q˙j(τ-)∫01Wj(x)Wk(x)dx+P∫01Wk(x)δd(x-xi)dx.
Making use of (10), the above equation reduces to:
(25)q˙j(τ+)=q˙j(τ-)+PWj(xi),
we now substitute (25) into (18):
(26)∑j=1N‍Wj(xi)[q˙j(τ-)+PWj(xi)]=-R∑j=1N‍Wj(xi)q˙j(τ-),
which can be simplified to obtain the following expression for P:
(27)P=-(1+R)∑j=1N‍Wj(xi)q˙j(τ-)∑j=1N‍Wj(xi)2.
The post-impact modal velocities can be obtained by substituting (27) into (25), which results in the following expression:
(28)q˙j(τ+)=q˙j(τ-)-(1+R)w˙(xi,τ-)∑j=1N‍Wj(xi)2Wj(xi),j=1,2,…,N.


### 2.2. Simulating Sticking Motions

During a chatter sequence the beam repeatedly impacts the obstacle and the time between two successive impacts continuously drops. The pre-impact velocity also continually drops and approaches zero.To switch to the sticking motion phase the time between two successive impacts is monitored and when this time falls below a certain threshold, the beam is assumed to stick with the obstacle. When the reaction between the obstacle and the beam vanishes the sticking motion is assumed to end.Figure2 illustrates the method of simulating a chatter sequence using a CoR method with Lagrange multiplier. In the figure, Δτi represents the time between successive impacts. When this time duration falls below a threshold value, ϵ, the beam is assumed to stick to the obstacle. To exactly satisfy the condition of sticking, as shown in Figure 2, when a chatter sequence is recognized the beam is brought to rest at the point of impact by applying an impulse using CoR = 0. The post-impact modal position and velocities can now be obtained from (17) and (28) using CoR = 0. The sticking phase is now simulated by applying a reaction force at the point of impact such that the sticking constraint is satisfied. The sticking phase ends when the reaction force changes its sign. We now derive the expression for this reaction force.Figure 2
Simulating a typical chatter sequence using CoR method with Lagrange multiplier.Letλ(τ) be the reaction force at the location of the impact during sticking. The modal equation, (13), is modified as follows by the addition of the reaction force:
(29)q¨j(τ)+2ζjωjq˙j+ωj2qj=αj(xf)Fsin⁡(ωτ)+Wj(xi)λ(τ),j=1,2,…,N.
The constraint equation resulting from the sticking condition is given by:
(30)w(xi,t)=∑j=1N‍Wj(xi)qj(τ)=d,
we differentiate (30) twice with respect to time:
(31)∑j=1N‍Wj(xi)q¨j(τ)=0,
substituting (29) into (31) we get:
(32)∑j=1N‍Wj(xi)(-ωj2qj(τ)-2ζjωjq˙j+αjFsin⁡(ωτ)+Wj(xi)λ(t))=0.
Solving the above equation for λ(t) we get
(33)λ(t)=-∑j=1N‍Wj(xi)αjFsin⁡(ωτ)∑j=1N‍Wj(xi)2-∑j=1N‍Wj(xi)(ωj2qj(τ)+2ζjωjq˙j)∑j=1N‍Wj(xi)2,
where λ(t) is the dynamic reaction at the location of the obstacle during sticking, which enforces the position constraint (30). Once the reaction approaches zero, the sticking motion ends and the modal equations of motion will be integrated with λ(t)=0. The process of simulating the vibro-impacting motion of a beam with sticking motions using Lagrange multiplier can now be summarized as below.(1)
Use the modal equations, (13), to integrate the states, {q(τ),q˙(τ)}, forward in time.(2)
Check for impact: if impact is detected use (28) to model the impact.(3)
Check for sticking: if time between two successive impacts is less than threshold,ϵ, the sticking phase begins, go to step 4. Otherwise continue with step 1.(4)
Integrate the modal equations, (32) which include the Lagrange multiplier (which represent the reaction at the pin) in the modal equations. The Lagrange multiplier, λ(τ), can be calculated using (33).(5)
Check for contact separation by monitoringλ(τ): if λ(τ) crosses zero sticking phase ends go to step 1; otherwise, continue with step 4.Before presenting numerical simulations we summarize the developments in this paper. In Section2 a mathematical model based on modal formulation has been developed for simulating nonimpacting motions of the beam. A modal form of CoR rule for simulating impacts was discussed in Section 2.1. The simulation of sticking motions using Lagrange multiplier is proposed in Section 2.2 within the framework of CoR rule.


## 2.1. Impact Modeling

Equation (13) when integrated with appropriate initial conditions can simulate the beam motion until the impact. Let τ- be the time just before the impact (pre-impact time). Pre-impact states, w(x,τ-) and w˙(x,τ-), are known. We have to find the states immediately after impact (post-impact), w(x,τ+) and w˙(x,τ+). Using the classical Newton CoR assumption for rigid impacts, the configuration of the beam is assumed not to change during impact [1] (as rigid impacts are assumed to occur in an infinitesimal time), we have
(16)w(x,τ+)=w(x,τ-).
The post-impact modal displacement initial conditions are
(17)qj(τ+)=qj(τ-),j=1,2,…,N.
We also know that the beam velocity at the location of impact is altered according to the classical CoR rule:
(18)w˙(xi,τ+)=-Rw˙(xi,τ-),
where R is the coefficient of restitution. We must now obtain the post-impact velocity distribution, w˙(x,τ+), satisfying (18). The procedure for obtaining the post-impact velocity distribution is explained in detail in Vyasarayani et al. [18]; we present a brief explanation here. Let us assume that at the instant of impact (τ-) an impulse of magnitude P acts on the beam at the impact location. The equation of motion of the beam can be written as:
(19)∂4w∂x4+c∂w∂t+∂2w∂τ2=Pδd(x-xi)δd(τ-τ-).
We assume the same pinned-pinned boundary conditions for the above PDE. Substituting the series solution given by (7) into the above PDE, we have
(20)∑j=1N‍d4Wj(x)dx4qj(τ)+∑j=1N‍Wj(x)q˙j(τ)+∑j=1N‍Wj(x)q¨j(τ)=Pδd(x-xi)δd(τ-τ-).
Integrating the above equation between τ- and τ+, we get
(21)∑j=1N‍d4Wj(x)dx4∫τ-τ+qj(τ)dτ+∑j=1N‍Wj(x)∫τ-τ+q˙j(τ)dτ+∑j=1N‍Wj(x)∫τ-τ+q¨j(τ)dτ=Pδd(x-xi)∫τ-τ+δd(τ-τ-)dτ.
The post-impact displacement distribution will remain the same as the pre-impact distribution, as the impact is assumed to occur in an infinitesimal amount of time, that is, w(x,τ+)=w(x,τ-). Hence we have qj(τ+)=qj(τ-) and thus the first integral in (21) becomes zero. Integrating the second term with respect to time, and from the definition of Dirac’s delta function, ∫τ-τ+δd(τ-τ-)dτ=1, (21) becomes
(22)∑j=1NWj(x)(q˙j(τ+)-q˙j(τ-))dt=Pδd(x-xi).
Now substituting for w˙(x,τ)=∑j=1N‍Wj(x)q˙j(τ), we get:
(23)w˙(x,τ+)=w˙(x,τ-)+Pδd(x-xi).
From the above equation we can say that the effect of applying an impulse on the Euler-Bernoulli beam results in a sudden change in velocity of the beam, immediately after the application of impact. The resulting velocity distribution after impact is nothing but a Dirac’s function in space added to the pre-impact velocity distribution. The scaling parameter P (magnitude of impulse) in (23) is to be found such that the post-impact velocity at the impact location is as described by (18). Multiplying both sides of (23) by Wk(x) and integrating over the domain results in:
(24)∑j=1N‍q˙j(τ+)∫01Wj(x)Wk(x)dx=∑j=1N‍q˙j(τ-)∫01Wj(x)Wk(x)dx+P∫01Wk(x)δd(x-xi)dx.
Making use of (10), the above equation reduces to:
(25)q˙j(τ+)=q˙j(τ-)+PWj(xi),
we now substitute (25) into (18):
(26)∑j=1N‍Wj(xi)[q˙j(τ-)+PWj(xi)]=-R∑j=1N‍Wj(xi)q˙j(τ-),
which can be simplified to obtain the following expression for P:
(27)P=-(1+R)∑j=1N‍Wj(xi)q˙j(τ-)∑j=1N‍Wj(xi)2.
The post-impact modal velocities can be obtained by substituting (27) into (25), which results in the following expression:
(28)q˙j(τ+)=q˙j(τ-)-(1+R)w˙(xi,τ-)∑j=1N‍Wj(xi)2Wj(xi),j=1,2,…,N.


## 2.2. Simulating Sticking Motions

During a chatter sequence the beam repeatedly impacts the obstacle and the time between two successive impacts continuously drops. The pre-impact velocity also continually drops and approaches zero.To switch to the sticking motion phase the time between two successive impacts is monitored and when this time falls below a certain threshold, the beam is assumed to stick with the obstacle. When the reaction between the obstacle and the beam vanishes the sticking motion is assumed to end.Figure2 illustrates the method of simulating a chatter sequence using a CoR method with Lagrange multiplier. In the figure, Δτi represents the time between successive impacts. When this time duration falls below a threshold value, ϵ, the beam is assumed to stick to the obstacle. To exactly satisfy the condition of sticking, as shown in Figure 2, when a chatter sequence is recognized the beam is brought to rest at the point of impact by applying an impulse using CoR = 0. The post-impact modal position and velocities can now be obtained from (17) and (28) using CoR = 0. The sticking phase is now simulated by applying a reaction force at the point of impact such that the sticking constraint is satisfied. The sticking phase ends when the reaction force changes its sign. We now derive the expression for this reaction force.Figure 2
Simulating a typical chatter sequence using CoR method with Lagrange multiplier.Letλ(τ) be the reaction force at the location of the impact during sticking. The modal equation, (13), is modified as follows by the addition of the reaction force:
(29)q¨j(τ)+2ζjωjq˙j+ωj2qj=αj(xf)Fsin⁡(ωτ)+Wj(xi)λ(τ),j=1,2,…,N.
The constraint equation resulting from the sticking condition is given by:
(30)w(xi,t)=∑j=1N‍Wj(xi)qj(τ)=d,
we differentiate (30) twice with respect to time:
(31)∑j=1N‍Wj(xi)q¨j(τ)=0,
substituting (29) into (31) we get:
(32)∑j=1N‍Wj(xi)(-ωj2qj(τ)-2ζjωjq˙j+αjFsin⁡(ωτ)+Wj(xi)λ(t))=0.
Solving the above equation for λ(t) we get
(33)λ(t)=-∑j=1N‍Wj(xi)αjFsin⁡(ωτ)∑j=1N‍Wj(xi)2-∑j=1N‍Wj(xi)(ωj2qj(τ)+2ζjωjq˙j)∑j=1N‍Wj(xi)2,
where λ(t) is the dynamic reaction at the location of the obstacle during sticking, which enforces the position constraint (30). Once the reaction approaches zero, the sticking motion ends and the modal equations of motion will be integrated with λ(t)=0. The process of simulating the vibro-impacting motion of a beam with sticking motions using Lagrange multiplier can now be summarized as below.(1)
Use the modal equations, (13), to integrate the states, {q(τ),q˙(τ)}, forward in time.(2)
Check for impact: if impact is detected use (28) to model the impact.(3)
Check for sticking: if time between two successive impacts is less than threshold,ϵ, the sticking phase begins, go to step 4. Otherwise continue with step 1.(4)
Integrate the modal equations, (32) which include the Lagrange multiplier (which represent the reaction at the pin) in the modal equations. The Lagrange multiplier, λ(τ), can be calculated using (33).(5)
Check for contact separation by monitoringλ(τ): if λ(τ) crosses zero sticking phase ends go to step 1; otherwise, continue with step 4.Before presenting numerical simulations we summarize the developments in this paper. In Section2 a mathematical model based on modal formulation has been developed for simulating nonimpacting motions of the beam. A modal form of CoR rule for simulating impacts was discussed in Section 2.1. The simulation of sticking motions using Lagrange multiplier is proposed in Section 2.2 within the framework of CoR rule.


## 3. Numerical Simulation and Discussion

A numerical code was developed in Matlab for simulating a pinned-pinned vibro-impacting beam (four modes were considered in the simulation). Matlab’s built-in integrator (ode15s) for solving ordinary differential equations along with the built-in event detection algorithm for detecting impacts was employed. We used an absolute and relative tolerance of1e-6 in the numerical simulations. For simulating the sticking motions the proposed Lagrange multiplier method within the framework of CoR rule was used. A threshold time for detecting chatter was chosen to be ϵ=1e-6 (later we show how this value can be selected). To validate our numerical code we first studied the free vibration response of the beam with CoR = 1. The initial conditions for the beam were chosen to be w(x,0)=3sin(πx), which corresponds to the first mode shape and w˙(x,0)=0. The beam is allowed to vibrate against an obstacle located at xi=0.4. The gap of d=0 was selected. The total energy of the beam is shown in Figure 3(a) and is conserved during the entire simulation (the beam impacts the obstacle 40 times during the simulation time), as required for perfectly elastic impacts. Figure 3(b) shows the phase plot of the beam and shows that the constraint (w(0.4,τ±)-d=0) is exactly satisfied at every impact.Free vibration response of the beam withw(x,0)=3sin(πx), R=1, and xi=0.4: (a) total energy of the beam, and (b) phase plot at the location of impact.
(a)
(b)


### 3.1. Choosing the Valueϵ

Figure4(a) shows a typical chatter sequence of the beam under external distributed harmonic forcing. Multiple impacts with reducing amplitudes are seen in Figure 4(b); the motion is similar to that of a bouncing ball falling under gravity on a perfectly rigid ground with R<1, where an infinite number of impacts occur before the ball comes to rest. The logarithm of time difference between the successive impacts (logeΔτ) corresponding to Figure 4(b) is shown in Figure 5. It can be seen that before the onset of sticking motions the time between successive impacts drops exponentially or linearly (almost) in the logarithmic scale. Now what if we stick the obstacle to the beam immediately when the chatter is detected (with ϵ=1e-6, the beam impacts the obstacle 19 times before it is made to stick); does it cause any significant change in the subsequent dynamics of the beam? To answer this question we have simulated the beam motion by choosing ϵ=1e-3 and ϵ=1e-4. These values are within the region of chatter sequence (see Figure 5). The phase plot of the beam motion at the location of impact for different values of ϵ is shown in Figure 6 and is evident from the figure that the overall motion of the beam (poststicking phase) remains the same irrespective of what epsilon value we choose, except the number of impacts the beam makes before it sticks to the obstacle are different for each ϵ. Once the chatter is detected (exponential decay of impact times between successive impacts) the beam can be made to stick to the obstacle without compromising on poststicking response behavior. The amount of computational time required to simulate the beam motion for different values of ϵ is shown in Figure 7. The comparison was made for 10 units of simulation time (the computational time was calculated by taking average computational time for 20 simulations). By using a lower value of ϵ=1×10-3 instead of ϵ=1×10-6, we were able to gain 16% of computational time. In order to characterize the nonlinear vibro-impacting continuous systems, several long-time integrations need to be performed while varying the detuning parameter, so any computational time gained will be beneficial while performing parametric studies.Predicted motion of the beam forF=76.8, ω=0.1ω1, R=0.7, d=0, ζj=0.05, xi=0.4, and ϵ=1×10-6: (a) beam displacement at the location of impact showing chatter and (b) close-up view near the end of chatter and the beginning of sticking phase.
(a)
(b)Figure 5
Impact duration between successive impacts for typical chatter sequence. The following parameters are used during the simulation:F=76.8, ω=0.1ω1, R=0.7, d=0, ζj=0.05, xi=0.4, and ϵ=1×10-6.Response of the beam for three different values ofϵ: (a) phase-plot at the impact location and (b) phase plot showing a close-up view near the impact.
(a)
(b)Figure 7
Comparison of computational time (for 10 units of simulation time) for different values ofϵ.


### 3.2. Regions of Sticking Motions

The sticking motions primarily occur due to wave phenomena in continuous systems. To emphasize the point consider this example; when the beam is subjected to impact at a location, the beam is abruptly brought to rest, while the rest of the beam continaes to move in the direction of the forcing. The impulse generated at the impact initiates a bending wave which will travel to the boundary and get reflected, thus causing further multiple impacts. Several parameters, like CoR, forcing amplitude, forcing frequency, and gap, can initiate sticking motions. As the intention of the paper is to demonstrate the capability of the developed theory for simulating vibro-impacting continuous systems, we will not discuss the specifics of motion, like studying the bifurcations and periodicity of motion, but rather we focus solely on sticking motions. Figure8 shows the Poincaré section at the location of impact with zero gap. In the frequency range of interest 0.5≤ω/ω1≤1.5, we found a large region where sticking motions (SMs) were present for the selected parameter values. Sticking motions are identified as zero rebound velocity impacts with finite time duration. It is evident from Figure 8 that the developed theory can be used to successfully simulate the vibro-impacting motion of the beam for a wide range of physical parameters.Figure 8
Poincaré section at the location of impact forF=72, R=0.7, d=0, ζj=0.05, xi=0.4, and ϵ=1×10-3 in the frequency range of 0.5≤ω/ω1≤1.5.


### 3.3. Control of Sticking Motions

The sticking motions can be controlled by changing different physical parameters, but we choose to vary the gapd for illustration. In this study, the same parameters that were used to generate Figure 8 are used, except we fix the frequency at ω=0.7ω1. Figure 9 shows the Poincaré section at the location of impact for different values of gap. As we increase the gap the sticking motions disappear around d≃0.6. Stable periodic double impact solutions (P-II), showing two impacts for every cycle of motion, were found over a wide range of gaps. Figure 10(a) shows the motion of the beam at the location of impact for a double impact motion for d=1.2. The corresponding phase space is shown in Figure 10(b). Other kinds of motions are also visible in Figure 9, like period-I, period-III, and also some bifurcations can be seen during sticking motions. A detailed parametric study of these motions in itself is an involved problem and is beyond the scope of the present work. The emphasis of this study was to show the capability of the developed theory to simulate a wide variety of impacting motions.Figure 9
Poincaré section at the location of impact forF=72, R=0.7, ζj=0.05, xi=0.4, and ϵ=1×10-3, ω=0.7ω1 and 0≤d≤1.8.Response of the beam forF=72, R=0.7, ω=0.7ω1, d=1.2,  ζj=0.05, xi=0.4, and ϵ=1×10-3: (a) displacement at impact location and (b) phase plot at impact location.
(a)
(b)


## 3.1. Choosing the Valueϵ

Figure4(a) shows a typical chatter sequence of the beam under external distributed harmonic forcing. Multiple impacts with reducing amplitudes are seen in Figure 4(b); the motion is similar to that of a bouncing ball falling under gravity on a perfectly rigid ground with R<1, where an infinite number of impacts occur before the ball comes to rest. The logarithm of time difference between the successive impacts (logeΔτ) corresponding to Figure 4(b) is shown in Figure 5. It can be seen that before the onset of sticking motions the time between successive impacts drops exponentially or linearly (almost) in the logarithmic scale. Now what if we stick the obstacle to the beam immediately when the chatter is detected (with ϵ=1e-6, the beam impacts the obstacle 19 times before it is made to stick); does it cause any significant change in the subsequent dynamics of the beam? To answer this question we have simulated the beam motion by choosing ϵ=1e-3 and ϵ=1e-4. These values are within the region of chatter sequence (see Figure 5). The phase plot of the beam motion at the location of impact for different values of ϵ is shown in Figure 6 and is evident from the figure that the overall motion of the beam (poststicking phase) remains the same irrespective of what epsilon value we choose, except the number of impacts the beam makes before it sticks to the obstacle are different for each ϵ. Once the chatter is detected (exponential decay of impact times between successive impacts) the beam can be made to stick to the obstacle without compromising on poststicking response behavior. The amount of computational time required to simulate the beam motion for different values of ϵ is shown in Figure 7. The comparison was made for 10 units of simulation time (the computational time was calculated by taking average computational time for 20 simulations). By using a lower value of ϵ=1×10-3 instead of ϵ=1×10-6, we were able to gain 16% of computational time. In order to characterize the nonlinear vibro-impacting continuous systems, several long-time integrations need to be performed while varying the detuning parameter, so any computational time gained will be beneficial while performing parametric studies.Predicted motion of the beam forF=76.8, ω=0.1ω1, R=0.7, d=0, ζj=0.05, xi=0.4, and ϵ=1×10-6: (a) beam displacement at the location of impact showing chatter and (b) close-up view near the end of chatter and the beginning of sticking phase.
(a)
(b)Figure 5
Impact duration between successive impacts for typical chatter sequence. The following parameters are used during the simulation:F=76.8, ω=0.1ω1, R=0.7, d=0, ζj=0.05, xi=0.4, and ϵ=1×10-6.Response of the beam for three different values ofϵ: (a) phase-plot at the impact location and (b) phase plot showing a close-up view near the impact.
(a)
(b)Figure 7
Comparison of computational time (for 10 units of simulation time) for different values ofϵ.


## 3.2. Regions of Sticking Motions

The sticking motions primarily occur due to wave phenomena in continuous systems. To emphasize the point consider this example; when the beam is subjected to impact at a location, the beam is abruptly brought to rest, while the rest of the beam continaes to move in the direction of the forcing. The impulse generated at the impact initiates a bending wave which will travel to the boundary and get reflected, thus causing further multiple impacts. Several parameters, like CoR, forcing amplitude, forcing frequency, and gap, can initiate sticking motions. As the intention of the paper is to demonstrate the capability of the developed theory for simulating vibro-impacting continuous systems, we will not discuss the specifics of motion, like studying the bifurcations and periodicity of motion, but rather we focus solely on sticking motions. Figure8 shows the Poincaré section at the location of impact with zero gap. In the frequency range of interest 0.5≤ω/ω1≤1.5, we found a large region where sticking motions (SMs) were present for the selected parameter values. Sticking motions are identified as zero rebound velocity impacts with finite time duration. It is evident from Figure 8 that the developed theory can be used to successfully simulate the vibro-impacting motion of the beam for a wide range of physical parameters.Figure 8
Poincaré section at the location of impact forF=72, R=0.7, d=0, ζj=0.05, xi=0.4, and ϵ=1×10-3 in the frequency range of 0.5≤ω/ω1≤1.5.


## 3.3. Control of Sticking Motions

The sticking motions can be controlled by changing different physical parameters, but we choose to vary the gapd for illustration. In this study, the same parameters that were used to generate Figure 8 are used, except we fix the frequency at ω=0.7ω1. Figure 9 shows the Poincaré section at the location of impact for different values of gap. As we increase the gap the sticking motions disappear around d≃0.6. Stable periodic double impact solutions (P-II), showing two impacts for every cycle of motion, were found over a wide range of gaps. Figure 10(a) shows the motion of the beam at the location of impact for a double impact motion for d=1.2. The corresponding phase space is shown in Figure 10(b). Other kinds of motions are also visible in Figure 9, like period-I, period-III, and also some bifurcations can be seen during sticking motions. A detailed parametric study of these motions in itself is an involved problem and is beyond the scope of the present work. The emphasis of this study was to show the capability of the developed theory to simulate a wide variety of impacting motions.Figure 9
Poincaré section at the location of impact forF=72, R=0.7, ζj=0.05, xi=0.4, and ϵ=1×10-3, ω=0.7ω1 and 0≤d≤1.8.Response of the beam forF=72, R=0.7, ω=0.7ω1, d=1.2,  ζj=0.05, xi=0.4, and ϵ=1×10-3: (a) displacement at impact location and (b) phase plot at impact location.
(a)
(b)


## 4. Conclusions

A new procedure for modeling sticking motions in vibro-impacting continuous systems has been proposed in the framework of the CoR method. The method does not require prior information of the mode shapes of the beam in sticking phase (when the beam is in contact with the obstacle) unlike mode-switching method. The sticking constraints are imposed exactly by calculating the closed-form expression for the constraint reactions using the Lagrange multiplier. A systematic method for selecting the sticking toleranceϵ based on time between successive impacts during a chatter sequence was proposed. A pinned-pinned vibro-impacting beam was numerically simulated using the developed theory to demonstrate the usefulness of the method.



---

*Source: 268595-2012-09-12.xml*