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.
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 () continuously decreases. Simulating the chatter sequence without assuming sticking when can be computationally prohibitive. For lumped parameter systems, Cusumano and Bai [23] proposed to monitor 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 Figure 1. The beam is subjected to a uniformly distributed harmonic force . The stop that limits the motion of the beam is located at . The governing partial differential equation (PDE) for the beam using Euler-Bernoulli beam theory, without considering the impact event, is as follows: where is the Young’s modulus, is the area moment of inertia, is the mass density, is the cross-sectional area, is the viscous damping coefficient per unit length, and is time. The pinned-pinned boundary conditions are enforced as follows: We scale the equation of motion by introducing the following parameters: The equation of motion in scaled form is as follows: The boundary conditions are The initial conditions are as follows: Using the expansion theorem, a solution to (4) is assumed to be of the following form: where and correspond to the th undamped mode shape and the corresponding modal coordinate for the beam, respectively. Substituting the above form of solution into (4), we get The above equation is multiplied by and then integrated over the spatial domain, which yields the following expression: We now evaluate the following integral: where is Kronecker’s delta, defined as: Since the amplitudes of the mode shapes are arbitrary, we choose the amplitude of mode shapes () to be . This procedure of choosing the amplitudes of mode shapes such that the right-hand side of (10) that becomes is called mass-normalization. By using the mass-normalized mode shapes, it can be shown that: 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 modes, are as follows: where and . 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: Multiplying both sides of (14) by , and then integrating over the domain, results in the following modal initial conditions:
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, and , are known. We have to find the states immediately after impact (post-impact), and . 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 The post-impact modal displacement initial conditions are We also know that the beam velocity at the location of impact is altered according to the classical CoR rule: where is the coefficient of restitution. We must now obtain the post-impact velocity distribution, , 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 acts on the beam at the impact location. The equation of motion of the beam can be written as: 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 Integrating the above equation between and , we get 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, . Hence we have 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, , (21) becomes Now substituting for , we get: 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 (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 and integrating over the domain results in: Making use of (10), the above equation reduces to: we now substitute (25) into (18): which can be simplified to obtain the following expression for : The post-impact modal velocities can be obtained by substituting (27) into (25), which results in the following expression:
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.
Figure 2 illustrates the method of simulating a chatter sequence using a CoR method with Lagrange multiplier. In the figure, 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.
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: The constraint equation resulting from the sticking condition is given by: we differentiate (30) twice with respect to time: substituting (29) into (31) we get: Solving the above equation for we get where 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 . 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, , 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 Section 2 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 of 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 (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 , which corresponds to the first mode shape and . The beam is allowed to vibrate against an obstacle located at . The gap of 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 is exactly satisfied at every impact.
(a)
(b)
3.1. Choosing the Value
Figure 4(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 , where an infinite number of impacts occur before the ball comes to rest. The logarithm of time difference between the successive impacts 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 , 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 and . 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 instead of , 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.
(a)
(b)
(a)
(b)
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. Figure 8 shows the Poincaré section at the location of impact with zero gap. In the frequency range of interest , 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.
3.3. Control of Sticking Motions
The sticking motions can be controlled by changing different physical parameters, but we choose to vary the gap for illustration. In this study, the same parameters that were used to generate Figure 8 are used, except we fix the frequency at . 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 . 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 . 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.
(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.
Disclosure
Some parts of this paper has appeared in “Effective Modelling of Sticking Motions in Vibro-Impacting Continuous Systems,” by C. P. Vyasarayani, S. S. Sandhu, and J. McPhee, Paper Number IMECE2010-39485, Volume 8: Dynamic Systems and Control, Parts A and B. Permission was obtained from ASME for reproducing the published material in this manuscript.
Acknowledgments
Funding of this research by the Natural Sciences and Engineering Research Council of Canada is gratefully acknowledged. Thanks to Dr. Kingshook Bhattacharyya of IIT Kharagpur for his technical inputs during the course of this study.