A walking roller chain

In this article, we study the motion of a vertically fixed rapidly spinning roller chain that, after being released, “walks” a certain distance along the floor. We construct a two-dimensional model of the roller chain. By numerically integrating its equations of motion, we make predictions on the distances travelled by the roller chain, the code of the simulation being available online in the supplements. Finally, we compare the predictions with experiments, discuss the results, and suggest topics for further research.


Introduction
The setup that we analyse in this article is as follows: A roller chain is placed around a vertically oriented sprocket, presumed to perfectly fit the roller chain. The circumference of the sprocket is smaller than the length of the roller chain, so a part of the chain hangs freely, with its bottom a small height above the ground. Initially, the chain and the sprocket are motionless. Next, we angularly accelerate the sprocket to a sufficiently high speed while keeping the roller chain on the sprocket (Fig. 1). After releasing the chain by lightly tapping it to the side, we notice that it traverses through a series of bounces a certain distance along the floor. This "walking" is the subject of our article. A video of the "walk" can be found online in the supplementary material.
In the video it can be seen that after being released, the roller chain falls a small height until the lower part of the chain collides with the floor. Due to the spinning of the chain, the impulse from the collision spreads across the chain and gives it an upward and forward velocity causing it to bounce. 1 The shape of the chain does not change significantly after the collision and the rotational speed of the chain is slightly reduced. The chain then moves freely until again the lower part collides with the floor, receives an impulse etc. This bouncing constitutes the roller chain's "walk", and it continues until the roller chain's angular speed becomes insufficient to sustain its vertical tilt, i.e. during its motion the roller chain's plane slowly tilts to the sides until finally it collapses on one side and stops from friction.
Previous studies of roller chains have predominantly investigated the mechanics of a chain within a chain drive system, its power transmission efficiency, longevity, durability and many other topics of practical significance [1][2][3][4]. Recently, successful attempts at modelling and simulating roller chains using multibody dynamics 2 have been developed [5][6][7], although still limited to chain drive motion. Currently there is no analysis of the motion of a free roller chain, rapidly spinning or otherwise, available in the literature. 3 The motivation of this article is to try to understand how this surprising and complex behaviour arises from the interplay of the roller chain with the environment, and to develop a model that captures this behaviour. Accordingly, in Section 2 we first develop a theoretical model of the roller chain. Then we present our numerical simulations and experimental measurements to test our model. We conclude with a detailed discussion of the agreement between theory and experiment.

Method
For our purposes, the most important thing to know about roller chains is that they are made of an alternating series of * e-mail: gpalle.phy@pmf.hr 1 By observing a slow-motion video of the collision between the roller chain and the floor it can be seen that the collision impulse creates two disturbances that spread along the chain in opposite directions in a wave-like manner. These disturbances appear in our numerical simulation as well, and it seems likely that they are the primary mechanism of the spread of the impulse. 2 Multibody dynamics is the study of the behaviour of a collection of rigid or flexible bodies interconnected in some arbitrary fashion by joints. 3 Many of the methods used in this paper, however, are well known in the real-time physics modelling and multibody dynamics literature, an introduction to which can be found in [8][9][10].
inner and outer chain links connected with plain bearings and that these plain bearings act as revolute joints that allow rotation around only one axis, restricting any other rotational or translational motion. Therefore, disregarding the slight twisting and bending of the joints, the roller chain may be adequately approximated as a planar object with revolute joints.

Roller chain model
In our model we take the inner and outer links to be ideal rigid bodies connected with revolute joints. If we circularly enumerate the chain links from 1 to n, taking odd indices to represent inner and even indices to represent outer links, then the value of the mass m i and moment of inertia I i of the ith chain link depends on the parity of i and is equal to if i is even: Since the roller chain's plane is close to the vertical during most of the "walk", we can justifiably fix it to the xyplane in our model, where x is the horizontal axis and y is the upward vertical axis (Fig. 2). Taking the gravitational field to be homogeneous with a gravitational acceleration g and using the coordinates of the centre of mass X i Y i and the orientation angle u i of the i th chain link, we can write the Lagrangian of the roller chain system as Revolute joints connect adjacent chain links, implying a set of time-independent holonomic constraints on the system. If ℓ denotes the distance between the joints of a chain link, we can write the holonomic constraints between the i th and (i + 1)th chain link as There are various ways of modelling the constraint forces that try to prevent the constraints from being violated. For numerical reasons, we treated the constraints by introducing very stiff springs between the joints. If k is the elastic constant of the springs, then the potential energy of the springs is Such a potential arises naturally if we think of the revolute joints as deformable bodies with a potential  minimum in the undeformed state. Understood in this way, an order of magnitude estimate of the elastic constant can be given by multiplying the length of a chain link with the Young's modulus of a typical metal: k ≈ 1 cm Â 100 GPa ¼ 10 9 N=m: Dissipative forces will be taken into account using generalized forces. If we introduce the generalized coordinates vector q ¼ ðX 1 ; . . . ; X n ; Y 1 ; . . . ; Y n ; u 1 ; . . . ; u n Þ T , we can write the corresponding generalized forces vector as with t i being the friction torque intrinsic to the joints and Q D j the generalized friction force arising from the deformation of the joints. The first 2n components of Q D j are the friction forces, while the last n components are the torques caused by these friction forces.
The rotational friction of the joints we take to be part viscous and part Coulomb's, with parameters b and t c , respectively (a common model for mechanical systems [11]). It is presumed that the lubrication and the pressure causing Coulomb's friction are the same across the roller chain. The discontinuity of the sign function in Coulomb's friction can significantly slowdown the numerical integration, because it disrupts the smoothness of the solution that is presumed in Taylor's theorem and used in the derivation of almost all numerical integration schemes. Thus we have smoothened the sign function with a hyperbolic tangent using a small v c parameter, the final expression being Without any damping, the stiff springs of the joints oscillate at high frequencies and cause unrealistic jerky motion of the chain model. Clearly, we have to account for damping during the deformation process of the joints, the particular model of damping being unimportant in our case, as long as it smoothens the motion. Even a very small damping was sufficient for smoothening the motion of the chain, while the motion itself depended unnoticeably on the amount of damping, most likely because after the transition to smooth motion, the deformation and the rate of change of the deformation became exceedingly small. The damping we took to be viscous with a parameter b and directed along the direction normal to the constraints: With q j as one of the X i , Y i or u i coordinates and Q j its respective dissipative generalized force component, the equations of motion can be written in component form as To write the equations of motion in vector form, we introduce the constraints vector C: . . . ; C x;n ; C y;1 ; . . . ; C y;n À Á T ; and using the vector of size n with all components equal to one 1, we introduce the gravitational acceleration vector g ¼ ð0; g1; 0Þ T of size 3n. We also define the 3n Â 3n massinertia matrix M as: Finally, we define the 2n Â 3n Jacobian matrix J ¼ ∂C=∂q of the constraints C with J ij = ∂ C i /∂ q j . The equations of motion can now be compactly written as For the numerical calculations it will be convenient to rewrite this second-order differential equation of 3n variables q as a first-order differential equation of 6n variables ðq; _ qÞ. Thus, the final equations of motion are d dt

Collisions
Collisions play a key role in the "walking" of the chain À they are the way the roller chain both gets its forward velocity and loses most of its energy. A comprehensive model of collisions must take into account both the tendency of the bodies to bounce off each other as well as the friction that tries to minimize slipping. One such model was proposed by Chatterjee and Ruina in reference [12]. The model is impulse-based, as it takes the collisions to be instantaneous, and algebraic, as it gives an algebraic relation between the pre-and post-collision quantities. Below we analyse a collision between two arbitrary two-dimensional bodies, that in the end of this subsection we specialize to those bodies appearing in the numerical simulation. Let A and B be two rigid objects that collide at a point P that is displaced a vector r A and r B from the centres of mass of objects A and B (Fig. 3), the positions related as and v B are their respective linear and angular velocities, then the velocities at point P of bodies A and B are and I B be the masses and moments of inertia of bodies A and B around their centres of mass and letn be the unit vector normal andt ¼ẑ Ân the unit vector tangential to the surface of body A at point P. If we mark the velocities after the collision with primes and take J ¼ J nn þ J tt to be the impulse transferred from body A to B, it follows that Using these equations we can get a relation between the transferred impulse and the change of the relative velocities of the bodies at P that, in general, has the following form: Choosingn andt as the basis for the matrix W and substituting equations (*) and (2) into the former equation, we get with the a coefficient defined as Now we specialize to the ChatterjeeÀRuina law, and define two auxiliary impulses: J 1 is the impulse necessary to bring the normal component of the relative velocity at point P to zero, while J 2 is the impulse necessary to bring the total relative velocity at point P to zero, using a frictional tangential component if necessary: In the following, m ∈ [0, 1] is the coefficient of friction, while e ∈ [0, 1] and e t ∈ [ À 1, 1] are the normal and tangential restitution coefficients, respectively. Using these auxiliary impulses we define the impulse J ðaÞ that annulates the relative tangential velocity and causes a bounce in the normal and tangential direction: In case this impulse violates Coulomb's law: |J t | mJ n , we project it on to the |J t | = mJ n cone using k to get J ðbÞ . It is easily seen that adding kðJ 2 À J 1 Þ to the impulse does not change the normal component of the relative velocity after the collision. Thus, the expression for J ðbÞ is J ðbÞ ¼ ð1 þ eÞJ 1 þ kðJ 2 À J 1 Þ; with k ¼ mð1 þ eÞJ 1;n jJ 2;t j À mðJ 2;n À J 1;n Þ : The ChatterjeeÀRuina collision law reads: The only collisions that appear in the numerical simulation are those between two chain links and between a chain link and an immovable object like the floor or sprocket. In the former case, we use the masses and moments of inertia of the corresponding chain links, while in the latter case we set m B = m i and I B = I i and take the limits m A ! ∞ and I A ! ∞ in all the preceding equations.

General remarks
Thus far we have developed a two-dimensional roller chain model, derived the equations of its motion (1) and determined a collision response 4 described by equations (2) and (3). Now we are in a position to discuss the "walking" of a roller chain within our theoretical model.
From the construction of the model it is apparent that the (theoretical model of a) chain as a whole accelerates downward because of gravity, that a quickly rotating chain has a tendency to continue rotating and sustain its shape because of the inertias of the individual chain links and that a collision of a chain link with the floor causes the chain as a whole to lose kinetic energy because of friction and bounce from the floor. The rotation of the chain, as well as the springs making the joints, aid in the transmission of the collision impulse to the whole chain. The internal frictions between the chain links cause an additional loss of kinetic energy. Consequently, our model seems to capture all the necessary aspects needed to manifest a "walking" motion that we see in reality.
Unfortunately, the complexity of the equations of motion make it difficult to be quantitative about any of these qualitative aspects of the motion of a roller chain. Only situations of high symmetry can be analysed analytically, e.g. the motion of a perfectly circular rotating chain in mid-air. In these situations, we can see some of the qualitative aspects of the roller chain's motion, but not all. What's more, these situations are too restrictive to make predictions for a comparison with an experiment. Hence, if we want to compare the theoretical model with an experiment, we will need to use numerical methods.

Numerical calculations
After trying various implicit and explicit methods in tackling the equations of motion (1), we have determined that high-order explicit methods with error estimators are to be preferred, as they allow for step-size adaptation. They were 40 times faster than implicit methods and 5 times faster than standard explicit methods, controlled for numerical error. 5 We used the eight-order RungeÀKuttaÀFehlberg 7(8) method that has a seventh-order error estimate, the Butcher tableau of which can be found in [14]. Our simulation also benefited from the odeint C++ library, available within the Boost C++ libraries package found in reference [15]. The code of the simulation is available in the supplementary material.
During the simulation we need to detect the collisions between the chain links, floor and sprocket, apply the appropriate impulses (Eq. (3)) and modify their linear and angular velocities (Eq. (2)) if they collide. For the broad phase collision detection 6 we used the sweep and prune algorithm, the details of which can be found in reference [16]. After finding a collision candidate, 7 we need to test whether and where the bodies intersect and determine the normaln and tangentialt unit vectors if they do. Given that this is simplest to do for circles, we model the chain links as two circles positioned at the joints 8 and the sprocket as a circle with a large friction coefficient and a small restitution coefficient. The floor we take to be the y = 0 plane.
The initial conditions in the simulation are obtained by placing the roller chain around the circular sprocket that is positioned at desired height above the ground, allowing the chain to fall on the sprocket and then slowly accelerating the sprocket to the desired rotational speed. Care is taken to ensure a steady rotating state before the release of the roller chain. The chain is released by instantaneously eliminating the sprocket.

Experimental method
The detailed behaviour of the spinning roller chain after being released (i.e. during its "walk") critically depends on the initial conditions, both in the experiment and simulation. Notwithstanding, we have found that the distance travelled by the roller chain varies only moderately with the length and initial rotational speed of the roller chain kept constant, therefore making the distance travelled our main quantitative measure for comparing the theory and experiment. A qualitative comparison is performed by comparing visualizations of the experiment and simulation and judging the physical soundness of the latter.
The experimental setup is as follows: A variablespeed drill is secured in place with a vice to a stand of variable height. A sprocket is attached to the drill, and the stand's height is set up so that the bottom of the roller chain almost touches the floor. With the roller chain placed on the sprocket, the drill is slowly accelerated to the desired rotational speed (Fig. 1) and released by gently tapping the spinning chain away from the drill with a metal rod. A barrier is placed between the sprocket and the drill to prevent the roller chain from slipping in the wrong direction and causing injury. After being released, the roller chain "walks" a certain distance Dx, defined to be the distance from the point directly below the centre of the sprocket to the estimated centre of mass of the resting roller chain. The "walks" whose trajectories significantly deviate from a straight line are disregarded as to make the comparison with the theoretical model more balanced. The experimentally varied parameters are the number of chain links n and the frequency of the sprocket n.
Various parameters are needed as input for the numerical calculations. The masses and dimensions of the chain links were measured directly, while the moments of inertia were calculated by making a model of its mass distribution (Fig. 4). The coefficients of friction between various parts of the system were either taken from a reference table or measured directly by uniformly dragging one body against the other with a dynamometer. The coefficients of restitution were looked up, as well as measured by dropping a chain link from a known height and evaluating the height of a vertical bounce using slowmotion footage. The rotational friction parameters b and t c were measured by constructing a pendulum around one of the joints, tracking its tip with a camera and fitting the parameters to the data. The rotational frequency n of the sprocket was determined using a slow-motion camera and the radius of the sprocket R was measured directly. The values of all the parameters used in the numerical simulation, together with their uncertainties, are given in Table 1. 5 An introduction to numerical integration of ordinary differential equations, as well as error control and step-size adaptation, can be found in reference [13]. 6 In an effort to reduce the complexity of collision detection from the Oðn 2 Þ complexity, in which you check every body against the other, broad phase algorithms have been developed. In them, space partitioning or sorting along an axis is used to discard the majority of pairs of bodies that need to be checked. 7 A collision candidate is any pair of bodies for which we suspect that they may be colliding. 8 In the majority of collisions, the middle non-circular part of a chain link can be disregarded (Fig. 4). A quantitative comparison of the simulated chain and the experiment is given in Figure 6. In the figure panels, Dx is either the measured or predicted distance travelled by the roller chain during its "walk" and n is the rotational frequency of the sprocket. Different panels show chains with varying numbers of chain links n, ranging from 60 to 180. Every experimental data point was obtained by averaging four measurements and calculating their standard deviation. Numerical points were calculated by simulating the "walk" for a given n and n 12 times and calculating the averages and standard deviations. For every simulation, a pseudorandom set of parameters was generated using Gaussian distributions with mean values and deviations as given in Table 1. 9 Except these parameters, in the simulation we have also varied the release time and height of the roller chain, thus resulting in further variations on the initial conditions.
In all the simulations marked red in Figure 6, the initial conditions were attained by simulating the acceleration of the sprocket and allowing the chain to reach a stationary state. In the n150 and n180 cases, this has resulted in an underestimate of the distance walked, and an inspection of the roller chain's motion showed that the cause of it was a premature collision of the chain with itself right after the second collision with the floor. In an effort to delay this collision, we have tried initializing the roller chain as spinning in a wider, more elliptical shape, resulting in the points drawn in green.
We have explored the sensitivity of the distance travelled on the numerical error tolerance used during integration. By increasing the precision thousand-fold, the distances travelled differed by up to 10%, given the same initial conditions. However, by taking a sample of 12 simulations, varied as described in the last paragraph, the averages differed by less than a few percent. As both low and high precision simulations have the same sensitivity to initial conditions and the same distributions of travelled distances, there is no physical reason to prefer one over the other within our (approximate) model. We have thus set numerical error tolerance to 10 À6 , rather than 10 À9 .
We have compared the qualitative behaviour of the roller chain, as seen on a slow-motion footage of the "walk", and a visualization of the numerical simulation (both available in the supplementary material) to conclude that our simulation captures all the qualitative aspects described in the introduction, except for the gradual lateral tilting of the roller chain's plain from the vertical. A graph depicting the dependence of the height and x component of the velocity of the centre of mass on time for a sample of simulations is given in Figure 5. The bouncing motion of the roller chain, abrupt changes in velocity due to collisions with the floor and energy loss caused by friction are all visible in Figure 5. For other chain link numbers and frequencies, the graph retains the same shape, with the number and size of the bounces varying.
In addition, using our model we have also investigated how the distance travelled depends on the various 9 The formula for the propagation of uncertainties based on linearization of the dependent quantity would have been inadequate in this case because the dependence of Dx on the various parameters appearing in the simulation is highly nonlinear in the range of one standard deviation, violating the assumptions of the formula.   Table 1. Increasing ℓ has resulted in a decrease of the distance travelled as tall chains topple more easily, while varying the masses and moments of inertia did not have an appreciable influence on the distance travelled. Neither did varying t c and b as any decrease in friction that resulted in less energy being dissipated came with a corresponding decrease in the rigidity of the roller chain, allowing it to lose its shape more easily. The coefficient of restitution of a floor-chain collision e fc we varied two orders of magnitude, and no dependence was seen. A decrease in the coefficient of friction of the floor m fc resulted in larger traversed distances because the chain was less likely to collide with itself right after a collision with the floor. The lower friction allowed longer sliding as well. Increasing the gravitational acceleration decreased the travelled distance, while the adjusted parameters in Table 1 had little effect on the distance travelled.

Discussion
From the results, it is evident that for the explored range of chain link numbers and frequencies, our model correctly predicts the distances travelled by the roller chain, given the proper initial conditions. The majority of simulated points fall within two standard deviations of the measured points. In light of the fact that there are no free parameters in our model to fit the data, that is all the parameters were measured independently, the agreement of theory with experiment is good. Our model is two-dimensional, and as such it never collapses on one side as seen in the experiment. Instead, the chain loses its shape in a more gradual manner, resulting in similar predictions on the distance travelled. An extension of our model in which the motion of the roller chain is decomposed to that of the centre of mass and that of the chain links within the chain's plane could take the gradual tilting of the plane into account. For this extended model, only six more degrees of freedom would be needed to represent the motion of the chain's centre of mass reference frame. The dynamics within the chain's plane would be the same as before except for the fictitious forces arising from the acceleration and rotation of the centre of mass reference frame.
Although we have varied the parameters that have the most pronounced influence on the distance "walked", the parameters of the floor and the roller chain used could also be varied. The only complication that arises with this variation is the difficulty in keeping all the other parameters constant when changing, for example, the mass of a chain link or the coefficient of friction of the floor.
An area of inquiry that we have not explored in any detail are the disturbances, mentioned in footnote 1, that appear when the roller chain collides with the floor. These wave-like disturbances can also be created by lightly tapping a roller chain that is rapidly spinning on a sprocket, as in Figure 1. Perhaps these disturbances could be explored with classical perturbation theory and it is possible that they influence the distance that the chain travels. Another phenomenon that consistently appeared in the experiment and has not been examined is the kink that develops at the bottom of the roller chain, also visible in Figure 1. To explain it, one might try searching for a steady-state solution of a roller chain whose motion is restricted by the sprocket.

Dead ends
We faced significant difficulties when attempting to use Lagrange multipliers to model the constraint forces, a detailed description of such a treatment can be found in the Constrained Dynamics lecture in [8]. In particular, the obtained equations of motion are stiff and as such require the use of implicit methods, which in turn need the calculation of a very complicated Jacobian. We have also attempted to reduce the number of generalized coordinates and the number of constraints. Such choices of generalized coordinates mostly lead to very complicated Lagrangians and equations of motion, while offering little computational benefit.

Conclusion
Using a simple model and computer simulations, we managed to replicate most qualitative aspects of the roller chain's motion such as the formation of the kink at the bottom while on the sprocket, the bouncing off the floor during the "walk" and the shockwaves created on impact with the floor. We also made quantitative predictions on the distances travelled during the "walk" that are in good agreement with experiment. We feel that this shows how simple ideas can go a long way in explaining challenging phenomena, keeping in mind their limitations.
We wish to thank Prof. Nikola Poljak for the discussions and aid offered during the writing of this article, as well as the referees for their constructive critique.