Supplemental Material for: “Dynamical elastic contact of a rope with the ground”

In Sec. I of this Supplement, we review the derivation of the dynamical equations of elastic rods and show of they simplify in the small-slope limit. In particular, we show how the hypothesis of a uniform tension closely follows from that approximation. We then express the governing equation in dimensionless form and derive, in Sec. II, the linear response of the rope to periodic actuation for arbitrary values of the bending stiffness. In Sec. III, we examine the effect of dissipation. Sec. IV provides details of the weakly nonlinear analysis. In Sec. V, we describe exact travelling wave solutions of the problem. These are used to implement the far-field radiation condition in the numerical algorithm, which we describe in Sec. VI.

In Sec. I of this Supplement, we review the derivation of the dynamical equations of elastic rods and show of they simplify in the small-slope limit. In particular, we show how the hypothesis of a uniform tension closely follows from that approximation. We then express the governing equation in dimensionless form and derive, in Sec. II, the linear response of the rope to periodic actuation for arbitrary values of the bending stiffness. In Sec. III, we examine the effect of dissipation. Sec. IV provides details of the weakly nonlinear analysis. In Sec. V, we describe exact travelling wave solutions of the problem. These are used to implement the far-field radiation condition in the numerical algorithm, which we describe in Sec. VI.

A. Nonlinear model
Mechanically, a rod designates a slender elastic body of cylindrical shape at rest. In the limit of a small aspect ratio, the extension of the rod can be neglected, to leading order, by comparison to transverse displacements (see the asymptotic treatment of the equations of nonlinear elasticity in [1] for a systematic justification). Denoting arc length along the rod bys and the local angle made by the tangent of the rod with the horizontal axis by θ(s), one has, for the coordinates (x(s, t), w(s, t)) of the center line of the rod, ∂x ∂s = cos θ(s), ∂w ∂s = sin θ(s).
Let T (s, t), N (s, t) and M (s, t) be the tension, shear force and torque at a section of coordinates in the rod. T and N are related to the horizontal and vertical forces as F x = T cos θ − N sin θ, F z = T sin θ + N cos θ.
Momentum conservation on an element ds of rod yields where ρ is the mass per unit length (ρ = ρ v A, where A is the transverse area and ρ v is the volume density). Equivalently, using Eq. (2), [2]: ρ cos θ ∂ 2 x ∂t 2 + sin θ ∂ 2 w ∂t 2 + g = ∂T ∂s − N ∂θ ∂s , ρ cos θ ∂ 2 w ∂t 2 + g − sin θ Next, conservation of angular momentum applied to the same element ds yields where I = (y 2 + z 2 )dydz is the geometrical inertia of the cross section of the rod. Furthermore, with E the Young modulus, one has the constitutive law [1] M = EI ∂θ ∂s . Hence, Eqs. (1), (4), and (7) form a complete differential system to determine x, w, T, N , and θ. * gkozyref@ulb.ac.be B. Small-θ limit The governing equation of the main paper is obtained in the limit θ 1. We first note in that limit that Eqs. (1) reduce to so that x ∼s and, hence, ∂ 2 x/∂t 2 ≈ 0. The first of Eqs. (4) then becomes implying that T is uniform along the rod and given by the value applied at the ends. Meanwhile, the remaining equations can easily be combined to give Except for the study of acoustic waves, the first term in the right hand side above is usually negligible and can safely be ignored. Hence, identifyings with x and defining the bending stiffness B = EI, one obtains In the above mathematical formulation, the hypothesis of a uniform tension T along the rod directly follows from the small-slope assumption, θ 1. When θ is not small, the horizontal acceleration of the element ds, the projection of gravity along the rod axis and the shear force all contribute to make T non uniform.
At x = 0, we have the boundary conditions w = Z 0 (1 + cos ωt) and B∂ 2 w/∂x 2 = 0, the later expressing [see Eq. (6)] that no torque is applied there. Next, at the contact point with the ground, x = x c (t), continuity requires that both w and ∂w/∂x vanish. In addition, since the ground exerts no torque, one also has B∂ 2 w/∂x 2 = 0.
It is convenient to rewrite Eq. (11) in dimensionless form for the following development. Let us adopt a unit length x 0 , a unit time t 0 , and a unit of altitude Z 0 . Hence, we write w(x, t) as where ξ and τ are the new, dimensionless, space and time. Eq. (11) then becomes Choosing t 0 and x 0 such that t 2 0 T /(ρx 2 0 ) = 1 and that gt 2 0 /Z 0 = 2, we eventually get Meanwhile, t 0 = 2Z 0 /g, x 0 = 2Z 0 T /ρg and the boundary conditions become, with Ω = ω 2Z 0 /g,

II. LINEAR ANALYSIS WITH FINITE BENDING STIFFNESS
Let us substitute in Eq. (14) the expansions We thus obtain Meanwhile, the boundary conditions become, at ξ = 0, On the other hand, to evaluate the boundary condition at ξ = ξ c , we Taylor-expand the functions near ξ c,0 to obtain Collecting all O(1) terms, we have to solve It is easy to check that solves both the differential equation for W 0 and the last four boundary conditions in (23). The remaining condition, W 0 (0) = 1, determines ξ c,0 : In the limit of a small bending stiffness, β 1, one easily finds from above that ξ c,0 ∼ 1 − β 1/2 + . . .. Moreover, the hyperbolic functions appearing in Eq. (24) only contribute significantly to W 0 (ξ) in boundary layers at the edges of the domain.
At O( ), the problem to be solved is with boundary conditions and To solve it, let us write yielding The general solution isŴ where µ and λ are given by The boundary conditionsŴ which determines the unknown constants A to D. These diverge if the determinant of the above matrix vanishes, that is if The first root of Eq. (34), given Eq. (25), is plotted as a function of β in Fig. 1. Let us look at two limiting cases: 1. β → 0: One has µ → Ω, λ → ∞, and ξ c,0 → 1. The resonance condition becomes sin Ω = 0, with Ω = π as the first resonance.
with Ω ≈ 4.45 as the first root, representing an increase of slightly less than 42% compared to π.
If the determinant of Eq. (33) does not vanish, one can solve it for A, B, C, and D and obtain W 1 (ξ, τ ). The displacement of the contact point ξ c,1 (τ ) then directly follows from Eq. (28).

III. LINEAR ANALYSIS WITH DISSIPATION
Let us examine the effect of dissipation on the linear response of the rope. For the sake of simplicity, we assume β → 0. Dissipation is introduced in the problem by modifying the equation of motion as Linear spectral response of the rope in the weak bending limit for an internal dissipation coefficient Γ = 0.3.
Above, Γ is a dimensionless internal friction coefficient. As β → 0 the remaining boundary conditions are As before, we expand the solution in powers of : At O( ), we get together with the boundary condition (21): The solution of Eqs. (38) is elementary: Hence, If the dissipation is weak (Γ Ω), then K ∼ Ω − iΓ and the general linear response to a driving function ∆Z(t)/Z 0 = A(ω) exp(iωt)dω becomes where γ = Γc/x 0 . Fig. 2 shows the modulus of the linear spectral response of the rope with finite Γ. Interestingly, large non-resonant frequencies can be more strongly amplified than even the first resonance. Hence, dynamical behaviors with high-frequency content can easily be excited.

IV. DETAILS OF THE WEAKLY NONLINEAR ANALYSIS
In this section, we provide details on the resolution of problems arising at higher orders in . As we have seen in Sec. II, the evaluation of the boundary conditions at the contact point is done by Taylor-expanding the solution near ξ = ξ c,0 . In order to avoid this, we introduce the shifted coordinate s = (x − x c (t))/x 0 = ξ − ξ c . Then, with one has, in the β → 0 limit, By construction, the contact point is now located exactly at s = 0. With β → 0, the boundary conditions are In order to evaluate the last of these conditions, we Taylor-expand W in the vicinity of s = −1: Since only one condition holds at s = −1 − d, only one Taylor expansion is required, thus simplifying subsequent manipulations. As before, we expand the variables in power of and collect terms of the same order in . We thus obtain an infinite sequence of problems, which we solve using Mathematica. We give the resolution details of the first orders. The leading order trivially yields The boundary condition at s = −1 suggests to write d 0 = D 0 cos Ωτ , with D 0 yet to be determined. Then the differential equation for W 1 becomes A particular solution associated with the right hand side above is Even though W 1,p vanishes at s = 0, its derivative does not. In order to satisfy ∂W 1 /∂s = 0 at s = 0, we must add to W 1,p a solution of the homogenous problem with the same temporal dependence. The general solution is W 1 = 2sD 0 cos Ωτ + (a 1 sin Ωs + b 1 cos Ωs) cos Ωτ and one quickly finds that a 1 = −2D 0 /Ω, b 1 = 0. Hence, in full agreement with Eqs. (40) and (41) when Γ = 0. The resolution at the following orders follows the same pattern, except that each time, new harmonics of the form cos jΩτ arise due to the nonlinear terms at O j . These harmonics are included in the function d i (τ ) and produce a set of inhomogeneous terms in the partial differential equation for W j (ξ, τ ). The particular solution are generally easy to find, and they are completed with homogeneous solution of the same temporal dependence in order to fullfil the boundary conditions at s = 0. Eventually, imposing the boundary condition at s = −1 determines the remaining unknown constants at this order. The calculation rapidly becomes heavy but can be managed with a symbolic software. Below, we summarize the results.
A. Near a resonance

Ω ≈ pπ
We now consider in detail the vicinity of the resonance Ω = pπ. Focusing on the limited expansion 3 j=0 j W j we let In this limit, we have Hence 3 j=0 j W j ∼ s 2 + (−) p pπs − sin pπs µ cos pπτ − 2 p 2 π 2 8µ 3 (2pπs − sin 2pπs) cos 2pπτ Keeping fixed, the last term above is the one that increases fastest as µ → 0. Indeed, if µ 1/2 , 3 j=0 j W j ∼ s 2 + 3 (−) p (pπ) 4 32µ 5 (pπs − sin pπs) cos pπτ + 1 3 (3pπs − sin 3pπs) cos 3pπτ , where the omitted terms are all smaller than the second one in the right hand side. After some manipulation, we obtain 3 j=0 j W j ∼ 1 p 2 π 2 p 2 π 2 s 2 + 3 (−) p (pπ) 6 32µ 5 (pπs − sin pπs) cos pπτ + 1 3 (3pπs − sin 3pπs) cos 3pπτ The function W (3) (S, φ) in Eq. (72) allows us to study the spatio-temporal expression for the rope near the contact point in the vicinity of all resonances Ω p = pπ. It is defined for S < 0, with S = 0 designating the contact point at all times. We now enquire when a new contact can be made with the ground for S < 0. By plotting the function W (3) (S, φ) as a function of its two arguments, one quickly concludes that a local spatial minimum can exist for S < 0 and that the value of that minimum is smallest when φ = 0. With this piece of information, we restrict our attention to W (3) (S, 0) and require The spatio-temporal dependence of W (3) (S, φ) is displayed in Fig. 3 for that value of k, illustrating grazing contact at φ = 0 and 2π.

Ω ≈ π/2
Let us now investigate the vicinity of Ω = π/2. To this end, we let At first sight, it would seem that it is enough to consider the limited expansion 2 j=0 j W j in order to obtain the conditions to produce a new contact, since the function W 2 in Eq. (56) diverges as 1/µ, through the term −Ω 2 8 sin 2 Ω cos 2Ω sin 2Ω (2Ωs − sin 2Ωs) cos 2Ωτ.
Following the same reasoning as in the previous section, we would find that, in the small-µ, with µ , However, the above expression is unable to properly predict touch down because such a function above can only touch down at s = −1, which is incompatible with the boundary condition that is imposed there. Hence, we must include higher order terms with higher spatial harmonics in order to ensure that a new contact takes place closer to s = 0.
Since W 3 remains bounded in the limit (78), we now consider the limited expansion where, as before, the omitted terms are all asymptotically smaller if µ is sufficiently small. Here, a new contact is made at s = −1/2 when cos 2πτ = −1 if 4 π 6 /(131072µ 3 ) ≈ 0.0796, that is if Note In the limit Ω−pπ = µ → 0, 4 W 4 actually dominates 3 W 3 if = O(µ 5/3 ) and, in all probability, higher-order terms grow even faster. This observation clearly indicates a deficiency of that limit. In the small-µ limit, one may conjecture from the first few orders of the calculation above that, at order n, n W n = O n µ −2n+1 . This suggests, as n becomes large, that the correct scaling near Ω = pπ is = O µ 2 , rather than the = O µ 5/3 announced in the paper. However, if one assumes this scaling at the start of the asymptotic analysis, one does not obtain a satisfactory solution indicating zero-velocity touch-down. This can be anticipated from the fact that all terms from the regular perturbation expansion become O(µ) and, hence are always too small compared to the static term s 2 . This is why, near Ω = pπ, we pragmatically truncate the expansion at third order. Even though the approach is not entirely satisfactory, the agreement with the numerically computed grazing bifurcation curves in the (ω, ) justifies it a posteriori.

V. TRAVELLING WAVES
Intuition suggests -and simulations and experiments confirm-that for sufficiently strong excitation, travelling waves are generated along the rope. Let us look for travelling wave solutions of Eq. (14). Writing Assuming that the longitudinal extent of the wave is 2 , we further have the boundary conditions Focusing on symmetrical waves, the general solution of Eq. (84) is The boundary conditions at ζ = are Defining they become a = −2β Finally, The first solutions of the transcendental Eq. (89) for C are {4.49341, 7.72525, 10.9041, 14.0662,. . . }. The travelling waves associated to these first 4 roots are represented in Fig. 4. The speed is indicating that small waves travel faster. That V depends on demonstrates the nonlinear character of the wave.

VI. NUMERICAL SIMULATIONS
A. Spatial discretization Space derivatives in Eq. (14) are discretized in the simplest possible way using second order finite difference. With a spatial discretization step ∆ξ, the state of the system is described by the time-dependent vector u(τ ) = {u i (τ )} where u i (τ ) = W (ξ i , τ ), ξ i = (i − 1)∆ξ. Eq. (14) is then approximated bÿ where overdots indicate time derivatives. Above, K is the differentiation matrix given by D − βD 2 , where D is the discrete laplacian: In order to impose the boundary conditions at ξ = 0, the first lines of Eq. (92) are replaced by In the right hand side of Eq. (92), q = {q i } is simply the gravitational force: q i = −2. Eq. (92) only holds for those values of ξ i for which the rope is not in contact with the ground -a situation that has to be re-evaluated at each time step.

B. The unilateral constraint
At each instant of time, there is a contact set that numerically defines the points of the rope that are considered in contact with the ground. δ is a small number, but not too small, in order to avoid spurious loss of contact due to discretization error (see later). If an element of the rope impacts the ground at some instant τ * , its velocity can undergo a discrete jump due to the impenetrability of the ground. Denoting the velocity just before and after the collision byu − andu + , respectively, one has the instant changeu Here λ is a vector of impulsive forces, non zero only at the indexes belonging to Σ τ and whose magnitude is not known a priori. The above abrupt change inu takes place in addition to the smooth evolution ofu given by Eq. (92), which takes into account the other forces that act simultaneously on each discrete segments of rope.

E. Influence of the contact detection threshold
Here, we discuss the effect of the parameter δ defined in Eq. (95) as the maximal altitude below which a point is deemed in contact with the ground. Fig. 6 illustrates the issue associated to setting its value. A small value, such as δ 0 may lead to an accurate evaluation of ξ c . However, because of the boundary condition of zero slope at contact, numerical noise can induce u i,k , together with neighboring discrete points of the rope, to parasitically cross that threshold. Such an event may result in large, abrupt variations of the numerical estimate of ξ c with an amplitude of several ∆ξ. Conversely, a larger value of δ such as δ 1 in Fig. 6, reduce the probability of parasitic excursion of ξ c , at the expense of accuracy. The determination of the numerical parameter δ is therefore a trade-off between accuracy and noise.
An illustration of this discussion is given in Fig. 7, which shows the time evolution of for three values of δ. Discontinuities of ξ c (τ ) signal the appearance of new contact points and, hence, wave emission. When wave generation starts, noise associated with the detection of ξ c becomes manifest, and is seen to increase with the lowering of δ. On the other hand, the largest value of δ in Fig. 7, while producing a smoother dynamics, underestimates the amplitude of oscillations of ξ c (τ ), as can be seen clearly at τ = 40. The latter is a general observation.
Given the noisy character of the dynamics in the wave radiating regime, particulary in association with the parameter δ, we now enquire whether the spikes in the bifurcation diagram in (Ω, ) can be a numerical artefact of the simulation algorithm. To this end we draw in Fig. 8 the bifurcation curves for different values of δ: 10 −2 , 10 −3 , 10 −4 , 10 −5 . We find that the computed bifurcation curves rapidly converge as δ decreases and that the spikes remain mostly at the same locations. In particular the dip at Ω ≈ 1.6 is well identifiable for all values of δ: This is the Ω = π/2 resonance discussed in Sec. IV A 2. Importantly, the curve obtained for the rather large value δ = 10 −2 , while vertically shifted compared to the other ones, still remains "spiky", with resonances mostly at the same frequencies as in the other curves. Hence, we may conclude that the spikes in Fig. 8 are not a numerical artefact associated to δ. Rather, the dynamical interpretation of these spikes as manifestation of nonlinear resonances is confirmed.
As discussed above, a large value of δ reduces parasitic excursions of ξ c and promotes a smooth dynamics, at the expense of accuracy, see Fig. 7. In our study, we performed most of our simulations with δ = 10 −3 .