Dynamics of ropes and chains: I. the fall of the folded chain

The problem of the paradoxical motion of the falling chain is considered. Laboratory and numerical experiments are performed with the initially folded configurations. The velocity and the acceleration of the falling tip are analysed. It is shown that if the acceleration of the chain tip exceeds the gravitational acceleration g, it results in the appearance of a sharp peak in the time evolution of the tip velocity. An approximate analytical formula describing the rising part of the peak is presented and reasons for its limited applicability are discussed.


Introduction
The dynamics of a rope is a purely classical problem. It proves to be relevant in explaining the motion of a number of physical objects such as: the whip, the fly line, the Indian rope, the space tether. Except for a few cases, such as the cracking whip whose motion is very fast, or the space tether, whose motion cannot be studied in a conventional laboratory, the dynamics of the rope is relatively easy to study experimentally. On the other hand, because of the complexity of the equations of motion, the theory of the problem is by no means simple. Not much can be said about the motion of the rope without applying numerical integration methods. As the algorithms of the numerical integration develop, the problem returns revealing its new, surprising aspects and details.
Most probably the first problem that has been considered both experimentally and theoretically is the motion of the whip, in particular, the cracking whip. Its long history has been recently reconstructed by Goriely and McMillen [1,2], who developed the contemporary theory of whip motion. Following the results of their bibliographic research, let us mention here the first experimental work done by Carriére [3] and the last experiment performed by Krehl et al [4].
Another context in which the dynamics of the rope proves to be essential is the casting of the fly line. The detailed description of the problem and its numerical analysis can be found in the latest papers by Gatti and Perkins [5,6], where one can find references to the earlier works by Spolek [7], Lingard [8] and Robson [9].
In both of the above-mentioned cases, one end of the rope is left free while the other end is fixed to a support moving along a desired trajectory. If the forced motion of the fixed end of the rope is periodic, one may expect the appearance of new phenomena such as the astonishing stabilization of the upside-down position of the rope. The latter phenomenon is known as the 'Indian rope trick'. One of the latest discussions on this can be found in [10].
A problem whose experimental study needs extraordinary technical means is the motion of the space tether. From the physical point of view, the space tether is a rope fixed by one of its ends to an orbiting space ship. Its dynamics are by no means simple, in view of the non-negligible gradient of the gravitational field and the non-inertial frame in which the motion takes place (see e.g. [11]).
Very often the physical model of the rope, convenient from the numerical point of view, is discrete. At zero elasticity and dissipation, such a model can also be seen as a model of the chain. Thus, in a natural manner, some problems of the rope dynamics can be seen as problems of the dynamics of chains. One of the most interesting problems of this kind is that of the paradoxical motion of the falling folded chain. Let us define it.
It has been known for a long time that the free end of a folded chain falls with an acceleration exceeding g. The falling chain problem was considered in detail by Schagerl et al [12]. The basic question asked in [12] was: why does the acceleration of the chain tip exceed g? The answer is by no means obvious, although some essential hints are provided by general considerations based on the energy conservation principle. As demonstrated by Schagerl et al, energy conservation leads to the appearance of an additional down-pulling force. In what follows, we reconsider the falling chain problem, asking somewhat different questions, using different initial conditions and comparing the experimentally observed motion with that simulated numerically. We check also the accuracy and the validity range of the approximate analytical formula, describing evolution of the chain-end velocity that can be derived using the energy conservation principle.

The discrete model of the chain and its equations of motion
The present paper is devoted to the falling chain and first we consider its equations of motion as a particular case of the equations of motion of a discrete model of the rope.
One can define a few discrete models of the chain; below we present one of them. Its equations of motion will be formulated for the case in which one of the ends is attached to the fixed support while the other one is free. Let us note that similar models have been considered before (see [12,13] and [10]).
Consider a chain moving in a gravitational field. There is only one force acting on it: the gravity force. Several assumptions will be made to simplify the model. First of all, we assume that the chain is constrained to move only in the two-dimensional plane; let it be the (x, y) plane. A chain of mass M and length L consists of n cylindrical rods (in what follows we shall refer to them as segments) with masses m i = m = M/n, i = 1, . . . , n, and lengths l i = l = L/n, i = 1, . . . , n. All the segments are rigid and cannot be deformed. Consecutive segments are connected by joints without elasticity and friction. Figure 1 shows the geometric representation of our model.
In order to introduce the equations of motion, the generalized coordinates must be specified; the coordinates must uniquely determine the state of the system. Our choice of the coordinates was motivated by our main requirement: an easy integrability of the resulting equations of motion. Following our predecessors, see e.g. [13], we decided to describe the system using angular coordinates indicating the inclination of the consecutive segments toward the x-axis.
The position of the first element is determined by the angle ϕ 1 . Similarly, the position of the second element is described by the angle ϕ 2 . The global configuration of the chain in the plane is uniquely expressed by ϕ i , i = 1, . . . , n. In the following discussion, the angles are referred to as generalized coordinates of the system. A generalized coordinate ϕ i indicates an angle between the ith element of a chain and the horizontal axis x.
The Cartesian coordinates of the ith mass centre (x i , y i ) can be written in terms of the generalized angular coordinates: Using the generalized coordinates, we shall derive the equations of motion. To start with we shall consider the energy of the system. The motion of the chain is considered as a combination of translational and rotational motions of its segments. Each segment has the moment of inertia I i = 1/12ml 2 , i = 1, . . . , n, calculated around the axis perpendicular to the (x, y) plane and passing through the centre of mass of the segment. Taking into consideration (1), the kinetic energy of the model is given by where() = ∂()/∂t. The potential energy of the ith segment is given by mgy i , where g is the gravitational acceleration. All in all, it may be expressed as The equations of motion of the falling chain are then governed by the system of Lagrange equations of the second kind: where L = T − U is the Lagrangian of the system [14]. Applying (2) and (3) to (4), we may find the equations describing the motion of a chain. Making all relevant simplifications, the motion of n connected rods are governed by a system of n differential equations of second order: where c i = cos(ϕ i ), c ij = cos(ϕ i − ϕ j ), s ij = sin(ϕ i − ϕ j ), q i = n − i + 1/2 and m i,j are elements of a symmetric and positive-defined mass matrix: For n = 2, we obtain the equations of motion of double rods: To make our considerations as general as possible, we find the dimensionless form of equation (5) by introducing the following dimensionless variable and parameter: where c denotes a characteristic velocity, e.g the velocity of sound in air. In the latter case, the results obtained can be easily used in the discussion of the cracking whip problem, where the velocity of sound plays an essential role. In the new dimensionless variables, the length of the chain and the characteristic velocity are equal to 1. Inserting (8) into (5) we derived the equation of motion in the non-dimensional form: where() from now denotes ∂()/∂τ. Thus, our goal has been achieved: we defined a discrete model of the chain and we derived its equations of motion.
Equations of motion formulated in this section cannot be solved analytically. Instead, their approximate solutions may be obtained using numerical methods. Introducing the new variables η i =φ i , i = 1, . . . , n, we rewrite the problem (9) as a first-order system of dimension 2n: The computation ofη i , i = 1, . . . , n, in the resulting system (10) requires the solution of a linear system. Then, the integration may be performed using a numerical integrator of Runge-Kutta type. We decided to use the RADAU5 algorithm designed for stiff problems and developed by Hairer and Wanner [15]. The method used in the code is based on a fifth-order Radau IIA scheme with an error estimator of order four.

Falling chain; initial configuration-broad catenary
The main aim of the present paper is to get a good insight into the dynamics of the falling chain. First of all, we would like to understand the details of its motion in the case when one of its ends is attached to the fixed support and the other end is allowed to fall freely. To reach the aim, we performed a series of simple laboratory and numerical experiments. Comparison of their results allowed us to obtain confidence in both the discrete model of the chain and its numerically integrated equations of motion.
In the laboratory experiments, we used a stainless steel chain consisting of 225 segments (see figure 2). Its length was equal 1 m and mass 0.022 kg. One of its ends was attached by a piece of thin thread to a firm support-a nail hammered into a wall. This made its rotation around the support point free. Let us denote the support point by O. The other end of the chain was attached with an easily removable fork-like element to another point of support. Let us denote it by P. Figure 3 shows the geometry of the experimental set-up.  Point P was located y = 0.75 m above and x = 0.22 m aside from the main support point O. The end attached to P was allowed to fall free only after all oscillations of the chain were damped out. The initial configuration that is formed as all the oscillations are damped out is close to the catenary curve. In what follows we shall refer to it as the broad catenary. A word of caution seems necessary here. From the rigorous point of view, the initial configuration formed in the experiments with chains built from a finite number of segments is not a catenary curve. One can see immediately that this is the case, considering the configuration of a chain built from an odd number of segments whose ends are fixed to points located at the same level and separated in the horizontal direction by the length of a single segment. In such a case the initial configuration consists of two exactly vertical pieces. This was the configuration used by Schagerl in his laboratory experiments [12]. We shall discuss it in the next section. We shall refer to it as the narrow catenary.
The initial velocity of all segments was zero. After removing the fork-like element, the chain was allowed to fall and its subsequent configurations were recorded with a digital camera. Figure 4 presents configurations of the falling chain found in the laboratory experiment. The brightness and contrast of the pictures recorded by the camera have been corrected to remove the unnecessary background and to visualize the best profile of the falling chain. The (x, y) coordination grid visible in the pictures was added to the pictures while processing.  between configurations found in the laboratory and numerical experiments is very good; this convinced us that the numerical simulations we perform describe properly the real dynamics of the falling chain. The discrepancy between conformations recorded in frames (e) will be discussed in the last section.
Performing the numerical simulations, we recorded all relevant data, in particular the trajectories in space of all segments of the simulated chain, their linear velocity and the potential and kinetic energies. In what follows we shall present and discuss their plots.
Let us start the discussion with the plot that shows the time dependence of the velocity at the end of the last segment, see figure 7. As one can see, it displays a sharp peak centred around the moment of time at which the falling chain becomes almost straight (see figures 5 and 6).
As indicated by Schagerl et al [12], the astonishing values reached by the velocity at the end of the falling chain can be qualitatively explained with an approximate analytical formula derived within a simplified model of the falling-chain dynamics. We assume in the model that the chain is divided into two linear parts: the motionless part attached to the support and the freely falling part. Initially, both parts are of equal length and they are motionless; the chain has a certain potential energy (one may assume that the chain has zero potential energy when it is freely hanging down) and zero kinetic energy. As the falling process starts, the potential energy diminishes and the kinetic energy grows. Obviously, in view of the energy conservation principle its sum remains a constant. The essential fact that makes further considerations different from the considerations carried out in the case of a compact falling body of a given mass, is that the mass of the falling part of the chain diminishes with time. Thus, the energy transferred from the potential to the kinetic form becomes squeezed into a smaller and smaller mass (the falling part of the chain). No wonder the acceleration of the chain end is greater than g. The energy-squeezing phenomenon, although not in such a pure form as in this simplified model, is present in reality. To visualize it, we plotted the numerically found distribution of the kinetic energy along the chain Such a definition of (s) makes its plot independent of the value of n. Looking at the plot one can easily see that the growing kinetic energy of the chain is indeed concentrated in this part of the chain, which is falling down; only a small amount of the kinetic energy is distributed along the part that is already hanging down. It seems to us that although the existence of the energy-squeezing phenomenon was expected, its clear, visual demonstration is new. The above qualitative reasoning carried out within the simplified model of the chain dynamics can be turned into a quantitative approximate calculation. In the calculation, the falling chain is treated as a variable mass system [16]. We assume the chain has a length L and that its mass M is distributed uniformly along it. One end of the chain O is fixed at y = 0. The other, falling end of the chain, is initially located at point P = (0, y 0 ) (see figure 3). We put x = 0, i.e. the free end of the chain is located exactly above the fixed end of the chain. To allow this, we assume that the number of segments from which the chain is assembled is infinite, i.e. the segments have infinitesimal length. In such conditions, one may assume that the falling chain is divided exactly into two vertical parts: the falling part (a) and the motionless, hanging downwards (b) (see figure 9). Dynamics of such a simple model can be solved analytically simply by applying the law of the conservation of energy.
The length and mass of the falling part of the chain are The length and mass of the motionless part are To find the potential energy of both parts, we must determine the vertical positions of their centres of mass: The most natural level in relation to which the potential energy of the whole chain can be determined is the level P 0 of its centre of the mass in the situation when the chain is freely hanging down. The vertical coordinate of this point is equal to −L/2. Thus, the potential energy of the falling part of the chain equals: Analogously, the potential energy of the motionless part of the chain is given by The total potential energy of the chain is now given by the formula: As can be easily seen, for y = −L, i.e. when the tip of the chain reaches its minimum position, U(−L) = 0. The kinetic energy of the falling part of the chain is given by Assuming that initially the free end of the chain is located at y 0 , we may formulate the law of conservation of energy for the falling chain: After straightforward simplifications we can find the formula describing the velocity v of the chain tip versus the actual position of its falling tip represented by y: v(y) = g y 2 0 + 2Ly 0 − y 2 − 2Ly L + y .
As seen from the formula, at the end of the falling process, when the whole chain becomes straight, y = −L and the velocity diverges.
To check the accuracy of such a simplified approach, we plotted the analytically predicted velocity v(y) together with the velocityv(y) found numerically in the simulations. See  (20),v(y) (red) is the velocity found in numerical simulations. Initial condition-broad catenary.
figures 10(a) and 10(b), where we plotted both the velocities and the relative difference between them defined as v(y) = |(v(y) −v(y))/v(y)|. Since the difference becomes significant only at the end of the falling process, we plotted the results for y < −0.5.
As seen in the figure, the differences are significant: (i) the maximum velocity of the realistic chain tip is finite and (ii) it appears earlier. As seen in figure 10(b) down to y = −0.9, the relative error of the simplified analytical prediction is smaller than 0.1 but for y < −0.9 the prediction fails. The main reason of that is rather clear: the idealized chain is built from the infinite number of segments, while the number of segments in a realistic chain is finite. Thus, at the end of the fall, when the falling part of the chain is reduced to a few segments, its discrete nature becomes essential and cannot be neglected. We shall discuss the problem in more detail below.

Falling chain; initial configuration-narrow catenary
As mentioned above, the shape of the initial configuration depends strongly on the initial position of the ends of the chain. When the ends are located at the same height and their horizontal separation is equal to the length of a single segment the configuration is particularly simple. For an odd number of segments the shape is exactly square: it consists of two vertically hanging pieces connected by a single segment; this was the initial configuration used in the experiments performed by Schagerl. Our type of chain does not allow such a configuration, but it can be easily defined in a numerical simulation, thus, we decided to repeat the Schagerl experiment simulating it numerically. Our aim was to check how the approximate analytical formula (20) fits its results. We assumed that the chain of length L = 1 m was built from n = 51 segments. In terms of the variables used in the analytical formula, the free end of the chain was initially located at y 0 = 0.
The shape of the initial configuration is shown in figure 11. As in the previous experiments, initial velocities of all segments were zero and at t = 0 one end of the chain was allowed to fall freely. Selected configurations of the falling chain found in the numerical simulations are presented in figure 12. As seen in the figure, the falling part of the chain remains straight almost to the last moment. On the other hand, the hanging part attached to the point of support develops visible oscillations. Consequently, a part of the growing kinetic energy of the chain is transmitted into it. As a result, the basic assumption of the analytical consideration presented above proves to be only approximate. Comparison of the chain-tip velocity found in the numerical experiments v(y) (red) and predicted by the approximate analytical theoryv(y) (black). Initial condition-narrow catenary. Figure 13 presents the time dependence of the chain-tip velocity found in the numerical experiment. The velocity can be compared with predictions of the approximate analytical formula (20). As seen in figure 14, the agreement is initially good, but below y = −0.8 significant differences appear. This happens since as the chain is falling a significant part of its energy becomes transmitted into the kinetic energy of the segments located in the hanging part of the chain. At the end of the falling process this happens to the last segment of the chain (see figure 15).
It seems worthwhile to point out another limit of the applicability of the approximate analytical formula (20). To analyse the problem better, let us consider first the fall of a ball attached to the support with a floppy, dissipation-free, massless and inextensible filament. The evolution of its velocity can be of course analysed using the energy conservation principle. As the ball reaches the lowest point, at which the filament becomes straight, the momentum of the ball changes its sign and the ball starts moving up. Obviously, the dynamics along this rising part of its trajectory can be also analysed using the energy conservation principle. What we find is that the flow of energy is simply reversed and the kinetic energy diminishes on account of the potential energy. The plot of the modulus of velocity versus time displays a symmetrical peak. The above considerations are somewhat trivial, but we recall them to confront their result with what happens in the falling-chain system. Formula (20) was derived using the energy conservation principle and, as we did in the falling ball case, one could be tempted to use it also after the moment of time, at which the falling chain becomes straight. Extending applicability of the formula after a moment one could conclude that having reached its lowest position, the end of the falling chain should start rising, reproducing eventually the initial configuration. This is certainly not observed in the numerical experiments. The rotational kinetic energy of segments cannot be easily transformed back into the translational energy leading to a smooth rising of the whole chain. What one observes instead is a very complex motion, where consecutive segments of the chain not only start moving up but also rotate. In general, the simulated motion is getting chaotic and, in view of the overlaps which appear, it loses the direct physical sense (see figures 12 and 13). Analysis of the problem requires a separate study.

Conclusions
It was known for a long time that the motion of the tip of the falling chain was not a motion with a constant acceleration. The question why this is the case has been posed by Schagerl. In [12] we find a convincing reasoning indicating that, in the particular configuration considered by the authors, the falling part of the chain is subject to an additional force pulling down the falling end and increasing its acceleration above the g value. Numerical simulations we performed revealed the exact profile of the velocity versus time dependence. The appearance of a sharp peak in the dependence can be seen as a result of the process in which almost all kinetic energy of the chain becomes squeezed into the shorter and shorter part of the chain. We have demonstrated that this is indeed the case. Observation of the phenomenon suggests the possibility of the derivation of an approximate analytical formula, whose accuracy, as we have shown, is rather limited. Moreover, after a moment of time, when the chain becomes straight, the formula fails completely.
Let us note that the existence of a very sharp peak in the time dependence of the velocity of the chain tip explains the origin of the discrepancy between the experimental and simulated configurations presented in frames (e) of figures 4 and 5. The procedure in the preparation of the figures was as follows. First, the configurations of the falling chain were recorded in laboratory experiments and the instants of time corresponding to the recorded configurations were determined. The values were subsequently used in the determination of the simulated configurations. Looking at figure 7 one can see that a tiny error in the determination of the time instant corresponding to the experimental configuration (e) may result in locating it at the wrong side of the velocity peak. The qualitative difference in the profiles of the experimental and simulated configurations visible in the frames (e) of figures 4 and 5 convinces us that we deal here with such a situation.
One of the most interesting questions one may ask while performing experiments with the falling chain is: at which length of the chain would the velocity of the tip exceed the velocity of sound? For practical reasons it is rather difficult to find a firm answer to the question, supported by a laboratory experiment. The answer is, on the other hand, easily provided by numerical simulations.
Consider a dissipation-free chain of length L built from n segments. Let v max be the maximum velocity reached by its tip. In view of the scaling properties of the equations of motion discussed above, changing the length of the chain toL > L while keeping the number of its segment fixed results in changing the maximum velocity of its tip to: Thus, demanding the maximum velocity to be equal to the velocity of the sound c = 340 m s −1 , we find the length L c of the chain which will be reached by this velocity: As calculations show, at the number of segments and the initial configuration discussed in section 3, the velocity of sound will be crossed when the length of the chain is larger than 21.5 m (v max ≈ 73.3 m s −1 ). For the number of segments and the initial configuration analysed in section 4, it will happen when the length is larger than 250 m (v max ≈ 21.5 m s −1 ). For a given length of a chain, the maximum velocity reached by the end of the chain strongly depends on the number of segments into which it is divided. This division has a direct physical sense in the case of a chain, but in the case of a continuous rope it is only a numerical tool and cannot influence results in such a dramatic manner. Numerical simulation of the motion of an elastic and dissipative rope is a challenging problem in which one can still find a lot of interesting details.