Impact of a Multiple Pendulum with a Non-Linear Contact Force

Abstract

This article presents a method to solve the impact of a kinematic chain in terms of a non-linear contact force. The nonlinear contact force has different expressions for elastic compression, elasto-plastic compression, and elastic restitution. Lagrange equations of motion are used to obtain the non-linear equations of motion with friction for the collision period. The kinetic energy during the impact is compared with the pre-impact kinetic energy. During the impact of a double pendulum the kinetic energy of the non-impacting link is increasing and the total kinetic energy of the impacting link is decreasing.

1. Introduction

Impact is a common and important phenomena in mechanical systems. The simplest impact analysis is based on the conservation of momentum and a kinematic coefficient of restitution defined by Newton. The kinematic coefficient of restitution is the ratio of the normal relative velocities after impact to the normal relative velocities before impact at the contact point [1]. Poisson divided impact into compression and restitution, and defined the kinetic coefficient of restitution as the ratio of the impulse during restitution to that during compression [2]. Kane and Levinson were able to extend the impact analysis and they showed for a double pendulum impacting a surface one can obtain energetically inconsistent results for friction using a kinematic coefficient of restitution [3]. Keller introduced the normal impulse as an independent variable which eliminated the normal force from the equations of motion [4]. Stronge defined a new energetic coefficient of restitution as the negative of the ratio of the work done by normal impulse during restitution to that during compression [5]. In contrast to the former two definitions, the energetic coefficient of restitution provides an energetically consistent solution for frictional collision. Stronge presented also a contact model using springs and the impulse at separation was obtained using the energetic coefficient of restitution [6]. All these previous research are using algebraic equations for the calculation of post-impact velocities. The dynamics of the systems with impacts are highly affected by the small changes of the coefficient of restitution which is estimated experimentally [7]. Impact responses of multibody systems have also been studied numerically in [8,9,10].

The double pendulum has been demonstrated to be chaotic for some configuration [11]. The expansion in separations of nearby trajectories was characterized by the Lyapunov exponent and an agreement between the results from experiment and numerical calculation was reached.

Feedback strategies have been introduced to control chaos of the impact of an inverted pendulum in [12]. The performance of the feedback was shown to be robust with respect to the variations of initial conditions and only minor information of the state of the system is required for such implementations. Litak et al. studied the effect of noise on the stability of the rotational motion of a parametric pendulum [13]. The use of the time-delayed feedback method in control of a parametric pendulum for wave energy harvesting has been explored in [14].

In practice, crane systems commonly suffer from the unexpected swings induced by the double-pendulum mechanism during transportation. In order to suppress the payload oscillations, the input shaping technique has been developed by smoothing the input commands generated by the human operator [15].

The pendulum system has been used in piezoelectric energy harvester (PEH) to amplify the energy harvested from human motion [16]. Experimental results showed that double pendulum based PEH dramatically improved the efficiency of the PEH, compared with PEH with cantilever beam and single pendulum. Experimental results were consistent with the numerical calculations which found that double pendulum system produced higher kinetic energy compared with the single pendulum. The performance of the conventional PID feedback controller and LQR optimal controller on double inverted pendulum system have been compared in [17]. Linearization and approximation were performed on the system equation using Jacobian matrix and an unstable equilibrium position. The simulation results showed that a desired response was obtained by implementing a LQR controller. Udwadia and Koganti proposed an analytical approach for controlling a multi-body inverted pendulum system without any linearazations or approximations [18]. An explicit expression of the control force was obtained by imposing a user-prescribed Lyapunov constraint on the system which minimized the control cost simultaneously. To compensate for the uncertain knowledge of the properties of the actual system, an additional control force was designed based on the generalization of the concept of a sliding surface. Numerical simulations verified that the system was able to move from an initial position to various inverted configurations using their methodology.

There has been a large number of publications that focused on the dynamics of double pendulum [11,12]. For the mechanical impact of a double pendulum previous researchers conducted studies to investigate the post-impact velocities. The developments have led to algebraic equations and three definitions of the coefficient of restitution (kinematic, kinetic, and energetic) [3,4,5]. Using a newtonian approach for some cases the total kinetic energy at the end of the impact is greater than the total kinetic energy at the beginning of the impact for a double pendulum [3]. Until now there is not a standard procedure for the calculations of the post-impact velocities of a double pendulum. We propose a differential model based on three phases of the impact: elastic compression, elasto-plastic compression, and restitution. For each impact intervals non-linear contact forces are developed. This model is different from previous impact algebraic and differential studies. The previous models did not consider the three periods and we propose a new elasto-plastic force and new permanent deformations for the impact of the double pendulum. With this method we solve the post-impact velocities without introducing a coefficient of restitution.

The energy loss of a mechanical system is a measure of utmost importance for impacting links and represents a new problem too. The kinetic energy during the impact is compared with the pre-impact kinetic energy. In this study we focus also on the kinetic energy of the impacting and non-impacting links. There are cases for which the kinetic energy of the non-impacting is increasing with respect to the initial impact kinetic energy. The total kinetic energy is decreasing during impact. This research provides a new strategy for analysis and design of impacting systems.

The reminder of the paper is organized as follows. In Section 2, the mathematical model of the general multiple pendulum in impact with a solid surface is formulated. The contact force expressions for the three impact phases are proposed in Section 3. The impact equations of motion are derived in Section 4. The numerical simulations and results for the impact of a double pendulum with a horizontal surface are discussed in Section 5.

2. Mathematical Model

Figure 1 represents a general multiple compound pendulum during the impact with a solid surface. The plane of motion will be designated the $\left(xy\right)$ plane. The y-axis is vertical, with the positive sense directed vertically upward. The x-axis is horizontal and is contained in the plane of motion. The z-axis is also horizontal and is perpendicular to the plane of motion. These axes define an inertial reference frame. The unit vectors for the inertial reference frame are ı, ȷ, and $\mathit{k}$. The links i are homogenous bars and have the lengths $L_i$ and the masses $m_i$. At O and $A_i$ there are pin joints. The mass center of links i is $C_i$. At a certain instant the free end B of the chain impacts the horizontal surface.

To characterize the configuration of the multiple pendulum, the generalized coordinates $q_i\left(t\right),i=1,\dots ,n$ are selected. The coordinate $q_i$ denotes the radian measure of the the angle between link i and the horizontal $x-$axis. The position vector and the velocity vector of the mass center $C_j$ of link j is

$$\begin{array}{ccc}\hfill {\mathbf{r}}_{C_j}& =& \left[\sum _{i=1}^j\left(L_{i-1}cos{q}_{i-1}\right)+\frac{1}{2}{L}_{j}cos{q}_j\right]\mathit{ı}+\left[\sum _{i=1}^j\left(L_{i-1}sin{q}_{i-1}\right)+\frac{1}{2}{L}_{j}sin{q}_j\right]\mathit{ȷ},\hfill \end{array}$$

(1)

where $L_0=0$. The velocity vector of the mass center $C_j$ of link j is

$$\begin{array}{ccc}\hfill {\mathbf{v}}_{C_j}& =& -\left[\sum _{i=1}^j\left(L_{i-1}{\dot{q}}_{i-1}sin{q}_{i-1}\right)+\frac{1}{2}{L}_j{\dot{q}}_{j}sin{q}_j\right]\mathit{ı}+\left[\sum _{i=1}^j\left(L_{i-1}{\dot{q}}_{i-1}cos{q}_{i-1}\right)+\frac{1}{2}{L}_j{\dot{q}}_{j}cos{q}_j\right]\mathit{ȷ}.\hfill \end{array}$$

(2)

The position vector and the velocity vector of the impact point B are

$$\begin{array}{ccc}\hfill {\mathbf{r}}_B& =& \sum _{i=1}^n{L}_{i}cos{q}_i\mathit{ı}+\sum _{i=1}^n{L}_{i}sin{q}_i\mathit{ȷ},\hfill \\ \hfill {\mathbf{v}}_B& =& -\sum _{i=1}^n{L}_i{\dot{q}}_{i}sin{q}_i\mathit{ı}+\sum _{i=1}^n{L}_i{\dot{q}}_{i}cos{q}_i\mathit{ȷ}.\hfill \end{array}$$

(3)

The kinetic energies of the link i is

$$\begin{array}{c}\hfill T_i=\frac{1}{2}{I}_{C_i}{\dot{q}}_i^2+\frac{1}{2}{m}_i{\mathbf{v}}_{C_i}·{\mathbf{v}}_{C_i},\end{array}$$

(4)

where $I_{C_i}$ is the mass moment of inertia of link i about the center of mass $C_i$. The generalized velocities are defined as

$$\begin{array}{c}\hfill u_i={\dot{q}}_i={\omega }_i,i=1,\dots ,n.\end{array}$$

(5)

The total kinetic energy of the multiple pendulum is

$$\begin{array}{ccc}\hfill T& =& T_1+T_2+\dots +T_i+\dots +T_n.\hfill \end{array}$$

(6)

3. Contact Force

The impact normal force at the contact point is calculated for three distinct periods: elastic phase, elasto-plastic phase, and restitution phase [19,20]. Figure 2 illustrates the normal contact force F during impact as a function of the normal elastic deformation $\delta$. The elastic compression phase starts from the beginning of the contact where $\delta =0$ and ends when the maximum elastic deformation, ${\delta }_e$, is reached. Next, the elasto-plastic compression phase begins and continues until the maximum deformation, ${\delta }_m$, is reached. At the maximum compression the relative velocity between the bodies is zero. The restitution phase starts at the maximum deformation and ends when there is no contact and a permanent deformation, ${\delta }_r$, is reached.

For the elastic phase the Hertz method is employed. The reduced modulus of elasticity, E, is calculated with

$$\begin{array}{c}E=\frac{E_r{E}_s}{E_s(1-{\nu }_r^2)+E_r(1-{\nu }_s^2)},\hfill \end{array}$$

(7)

where $E_r$ is the modulus of elasticity of the impacting rod, $E_s$ is the modulus of elasticity of the flat surface, ${\nu }_r$ is the Poisson ratio of the impacting rod, and ${\nu }_s$ are the Poisson ratio of the flat surface. The reduced radius, R, is calculated with

$$\begin{array}{c}{R}^{-1}=\frac{1}{R_r}+\frac{1}{R_s},\hfill \end{array}$$

(8)

where $R_r$ is the radius of curvature of tip of the impacting rod and $R_s=\infty$ are the radius of curvature of the surface.

I. For the elastic phase the normal impact force is

$$F=\frac{4E{R}^{0.5}{\delta }^{1.5}}{3},$$

(9)

where $\delta$ is the normal elastic deformation

$$\begin{array}{c}\hfill \delta ={\mathbf{r}}_B·\mathit{ȷ}-\sum _{i=1}^n\left(L_{i}sin{q}_{i0}\right).\end{array}$$

The initial impact position is given by the angles $q_{i0}$. The elastic phase is ending when $\delta >1.9{\delta }_c$. The critical deformation, ${\delta }_c$, is calculated with

$${\delta }_c=R{\left(\frac{\pi C_j{S}_y}{2E}\right)}^2,$$

where ${\displaystyle C_j=1.295e^{0.736\nu }}$ is a critical yield stress coefficient, and $S_y$ is the yield strength of the weaker material.

II. For the elasto-plastic phase the force is calculated with a modified Jackson and Green expression [19,20]

$$\begin{array}{c}\hfill {\displaystyle F=P_c\left\{{\displaystyle e^{-0.17{(\delta /{\delta }_c)}^{5/12}}{\left(\frac{\delta }{{\delta }_c}\right)}^{1.5}+\frac{4H}{C_{j}Sy}\left[1-e^{(-1/78){\left(\delta /{\delta }_c\right)}^{5/9}}{\left(\frac{\delta }{{\delta }_c}\right)}^{1.1}\right]}\right\},}\end{array}$$

(10)

where

$$\begin{array}{c}{\displaystyle B=0.14e^{23S_y/E},a=\sqrt{{\displaystyle R{\delta }_c{\left({\displaystyle \frac{\delta }{1.9{\delta }_c}}\right)}^B}},\frac{H}{S_y}=2.84-0.92\left[1-cos\left(\pi \frac{a}{R}\right)\right],}\\ {\displaystyle P_c=\frac{4}{3}{\left(\frac{R}{E}\right)}^2{\left(\frac{\pi C_j{S}_y}{2}\right)}^3.}\end{array}$$

III. For the restitution an elastic impact force is considered

$$\begin{array}{c}\hfill F=\frac{4E{R}^{0.5}{(\delta -{\delta }_r)}^{1.5}}{3}.\end{array}$$

(11)

The permanent deformation, ${\delta }_r$, is given by [19,20]

$$\begin{array}{c}{\delta }_r={\delta }_m\left\{0.8\left[{\displaystyle 1-{\left(\frac{{\delta }_m/{\delta }_c+5.5}{6.5}\right)}^{-2}}\right]\right\}.\end{array}$$

(12)

The maximum deformation at the end of the compression phase is ${\delta }_m$. The radius of curvature of the tip for the restitution phase is

$$\begin{array}{c}\hfill R_r=\frac{1}{{\left({\delta }_m-{\delta }_r\right)}^3}{\left(\frac{3P_m}{4E}\right)}^2,\end{array}$$

where $P_m=F(\delta ={\delta }_m)$ is the maximum elasto-plastic force.

4. Impact Equations

During the impact the Lagrange equations of motion are used

$$\begin{array}{c}\hfill \frac{d}{dt}\left(\frac{\partial T}{\partial {\dot{q}}_i}\right)-\frac{\partial T}{\partial q_i}=Q_i,i=1,\dots ,n\end{array}$$

(13)

where $Q_i$ are the generalized active forces. The friction force at the contact point B is

$$\begin{array}{c}\hfill {\mathbf{F}}_f=-\mu \left|\mathbf{F}\right|\frac{({\mathbf{v}}_B·\mathit{ı})\mathit{ı}}{|{\mathbf{v}}_B·\mathit{ı}|},\end{array}$$

(14)

where $\mu$ is the kinematic coefficient of friction, $\mathbf{F}=F\mathit{ȷ}$, F is the normal impact force defined in Equations (9)–(11) for elastic compression, elasto-plastic compression and restitution phase, respectively.

The gravitational forces on link i, at the mass center $C_i$ is

$${\mathbf{G}}_i=-m_{i}g\mathit{ȷ},$$

where g is the gravitational acceleration. The generalized active forces associated to $q_i$ are

$$\begin{array}{ccc}\hfill Q_i& =& {\mathbf{G}}_i·\frac{\partial {\mathbf{v}}_{C_i}}{\partial u_i}+(\mathbf{F}+{\mathbf{F}}_f)·\frac{\partial {\mathbf{v}}_B}{\partial u_i},i=1,2,\dots ,n.\hfill \end{array}$$

(15)

For the dynamics of the system before and after the impact equations given by Equation (13) are used with $\mathbf{F}=\mathbf{0}$.

Numerical methods for kinematic chain dynamics have been studied extensively in the engineering and mathematics. To simulate complex mechanical systems, such as robots and walking machines, numerical methods based on ODEs and/or DAEs are used. To solve the discontinuity related to friction a regularized function can be used [21,22]. For each impact phase of the mechanical system, the contact forces were generated and feed into a system of ordinary differential equations. Numerical MATLAB methods were used to get an accurate approximate solution to the differential equations. The numerical problems have unique solutions but the continuous problem might have multiple solutions. Repeated simulations with random disturbances should be performed and sets of possible solutions can be determined [21].

5. Application and Results

The method in this study can be applied to a system with n links. For the numerical application we selected the impact of a double pendulum with a horizontal surface. Figure 3 shows a double pendulum with homogeneous links made of steel with the length $L_1=L_2=L$, diameter $d=2R$ and mass $m_1=m_2=m$. The impact flat surface is also made of steel. The density of the material is $\rho$, the Young’s modulus is E, the Poisson ratio is $\nu$, and the yield strength is $S_Y$. The kinematic coefficient of friction is $\mu$. The initial impact angles of the links are $q_{10}$, $q_{20}$. The impact point is defined by the position vector ${\mathbf{r}}_{B0}=\left(L_{1}cos{q}_{10}+L_{2}cos{q}_{20}\right)\mathit{ı}+\left(L_{1}sin{q}_{10}+L_{2}sin{q}_{20}\right)\mathit{ȷ}.$ The initial impact angular velocities of the links are ${\omega }_{10}$, ${\omega }_{20}$. The material properties, geometries and initial conditions used for the simulation are shown in

Table 1. Material properties, geometries for the double pendulum and the impact flat surface.

Double Pendulum		Flat
$\rho$	7800 (kg/m3)	$\rho$	7800 (kg/m3)
E	210 (GPa)	E	210 (GPa)
$\nu$	0.29	$\nu$	0.29
$S_Y$	1.12 (GPa)	$S_Y$	1.12 (GPa)
$\mu$	0.2	$\mu$	0.2
L	1 (m)
R	0.005 (m)
m	1 (kg)

and

Table 2. Initial conditions of the double pendulum.

Link 1		Link 2
$q_{10}$	$-{70}^{\circ }$	$q_{20}$	$-{60}^{\circ }$
${\omega }_{10}$	$-0.1$ (rad/s)	${\omega }_{20}$	$-0.2$ (rad/s)

.

The analytical expressions for the Lagrange equations of motion are obtained using symbolic MATLAB. The numerical solution of the non-linear ordinary differential equations is calculated with the MATLAB ode45 function. The end of elastic compression, maximum compression, and the end of restitution are detected with MATLAB event functions.

Figure 4 shows the generalized coordinates $q_1$ and $q_2$ during impact. The magnitude of the angle $q_1$, of the link 1, is increasing, as seen in Figure 4a. The magnitude of the angle $q_2$, of the impacting link 2, is increasing initially and then is decreasing, as depicted in Figure 4b. The maximum angle of the impacting link is obtained during the elasto-plastic phase.

Figure 5 shows the angular velocities ${\omega }_1$ and ${\omega }_2$ during impact. The non-impacting link 1 has a pure rotational motion and is connected with the impacting link 2 at a frictionless revolute joint, A. The impact force does not act directly on link 1. The motion of link 1 is influenced by the joint reaction force at A (the reaction force of the impacting link 2 on link 1) and the weight $G_1$. The magnitude of the angle $q_1$ is increasing during the impact and the magnitude of the angular velocity of link 1 will increase during the impact [3,23]. Figure 5a shows the increase of the magnitude of the angular velocity of link 1, ${\omega }_1$.

The impacting link 2 has a general motion (translation and rotation) and the normal impact force at B, the tangential friction force at B, the weight $G_2$, and the joint reaction force of link 1 on link 2 at A are considered for this element. The angular velocity ${\omega }_2$ has a complex motion: due to the impact and friction force the angular velocity is decreasing until zero, then is changing the sign, and next is increasing [23]. For the given geometrical configuration, Figure 5b shows the angular velocity of the impacting link 2, ${\omega }_2$. The angular velocity of the impacting link is changing the sign during the elasto-plastic phase and its magnitude is increasing to a final value larger than the initial value. Simulations for another configuration of the double pendulum with $q_{10}=-{90}^{\circ }$, $q_{20}=-{78}^{\circ }$ shows that final angular velocity is ${\omega }_{2f}=0.178$ rad/s and it is less than the initial magnitude value (${\omega }_{20}=-0.2$ rad/s). For this case the angular velocity ${\omega }_2$ changes sign at maximum compression. The evolution of the angular velocity of the impacting end is dependent on the initial conditions.

The velocity of the contact point B, during impact, is represented in Figure 6. The normal velocity $v_{By}$ is decreasing and for the maximum compression $v_{By}=0$, as seen in Figure 6a. The tangential velocity of the impact point, $v_{Bx}$, is shown in Figure 6b. The tangential velocity is decreasing and changes the sign during restitution phase.

Figure 7a depicts the deformation $\delta$ during impact. The maximum deformation is at the end of compression. Figure 7b shows the normal force F during impact. The normal force is maximum at maximum compression.

Figure 8a shows the increase of the kinetic energy of the non-impacting link 1 during collision. The kinetic energy of link 1 is dependent on the square of the angular velocity ${\omega }_1^2$ and that is why the kinetic energy of link 1 will increase with the increase of the magnitude of ${\omega }_1$, as seen in Figure 5a and Figure 8a. For the impacting link the kinetic energy is decreasing as shown in Figure 8b. The total kinetic energy, shown in Figure 8c, is decreasing during impact.

Figure 9 shows the total kinetic energy for $q_{10}=-{70}^{\circ }$, $q_{20}=-{60}^{\circ }$, and ${\omega }_{10}=-0.1$ rad/s. The initial angular velocities of the impacting link are ${\omega }_{20}=-0.2,-0.3,-0.4,-0.5$ rad/s. The energy loss is defined as the difference between the initial kinetic energy and the final kinetic energy, $\Delta T=T_0-T_f$. The energy loss, $\Delta T$, is increasing with the magnitude of the initial angular velocity, ${\omega }_{20}$, of the impacting link: $\Delta T=0.005$ J for ${\omega }_{20}=-0.2$ rad/s, $\Delta T=0.010$ J for ${\omega }_{20}=-0.3$ rad/s, $\Delta T=0.016$ J for ${\omega }_{20}=-0.4$ rad/s, and $\Delta T=0.025$ J for ${\omega }_{20}=-0.5$ rad/s. For this case impact duration is decreasing with respect to the magnitude of the initial angular velocity, ${\omega }_{20}$. The energy loss is proportional with the permanent deformation ${\delta }_r$ as seen in Figure 10a: $\Delta T=0.005$ J for ${\delta }_r=5.49\left({10}^{-6}\right)$ m, $\Delta T=0.010$ J for ${\delta }_r=8.64\left({10}^{-6}\right)$ m, $\Delta T=0.016$ J for ${\delta }_r=1.25\left({10}^{-5}\right)$ m, and $\Delta T=0.025$ J for ${\delta }_r=1.55\left({10}^{-5}\right)$ m. Figure 10b shows the increase of the permanent deformation with the increase of the magnitude of the initial angular velocity, ${\omega }_{20}$. This is one of the reason for the increase of the energy loss with the magnitude of the initial angular velocity of the impacting link.

Figure 11 depicts the total kinetic energy of the double pendulum, during collision, for different initial impact angles of the impacting link: $q_{20}=-{50}^{\circ }$, $q_{20}=-{55}^{\circ }$, $q_{20}=-{60}^{\circ }$, and $q_{20}=-{65}^{\circ }$. Figure 12 shows the permanent deformation as a function of initial angle of the impacting link: ${\delta }_r=5.30\left({10}^{-6}\right)$ m and $\Delta T=0.005$ J for $q_{20}=-{50}^{\circ }$ m, ${\delta }_r=5.33\left({10}^{-6}\right)$ m and $\Delta T=0.005$ J for $q_{20}=-{55}^{\circ }$ m, ${\delta }_r=5.49\left({10}^{-6}\right)$ m and $\Delta T=0.005$ J for $q_{20}=-{60}^{\circ }$ m, and ${\delta }_r=5.82\left({10}^{-6}\right)$ m and $\Delta T=0.005$ J for $q_{20}=-{65}^{\circ }$ m. The permanent deformation is not influenced by the initial angle of the impacting link, and the energy loss is not influenced by the initial impact angles of the impacting link, $\Delta T\approx 0.005$ J. The duration of the impact is increasing with the magnitude of the initial impact angle.

Another reason of the increase of the kinetic energy loss with the magnitude of the impact angular velocity as seen in Figure 9 is that the kinetic energy is dependent on the square of the angular velocity, ${\omega }_2^2$. The kinetic energy is also dependent on the square of the velocity of the center of mass of link 2, $v_{C_2}^2$. The impact angle, $q_2$, has a small influence on the square of the velocity of the center of mass of link 2 and that is why the kinetic energy and the kinetic energy loss are less dependent on the impact angle as depicted in Figure 11.

6. Conclusions

This article presents a method to solve the impact of a kinematic chain in terms of a nonlinear contact force. The non-linear differential equations for the impact are obtained using Lagrange equations. The contact force is calculated for elastic compression, elasto-plastic compression, and elastic restitution. The final impact time is obtained from the permanent deformation of the material. The angular velocities of the double pendulum are increasing during the impact. The tangential linear velocity of the contact point is changing the sign during restitution phase. The kinetic energy of the non-impacting link is increasing during the impact. The total kinetic energy of the pendulum is decreasing during the impact period. To validate the proposed model, experimental investigation of the impact of the double pendulum will be developed. A high-speed camera can be used to capture the motion of the double pendulum before, during and after impact by placing markers on the links. The permanent deformation after impact can be scanned and measured with an optical profilometer. The presented approach can be applied to complex dynamical, systems such as robots, walking machines, and animal locomotion.