Propagation of elastic solitons in chains of pre-deformed beams

Following the seminal numerical experiment of Fermi–Pasta–Ulam–Tsingou [1], which was related by Zabusky and Kruskal to the propagation of solitons [2], a variety of model equations, solution methods and experimental platforms have been developed to investigate the dynamics of discrete and nonlinear one-dimensional mechanical systems across many scales [3–8]. At the macroscopic scale, propagation of solitary waves has been observed in a variety of nonlinear mechanical systems, including chains of elastic beads [9–14], tensegrity structures [15], origami chains [16], wrinkled and creased helicoids [17] and flexible architected solids [18–21]. Moreover, it has been found that even at the molecular scale solitons affect the properties of a variety of one-dimensional structures, including macromolecular crystals [22], polymer chains [23–26], DNA and protein molecules [27–30]. Since detailed experimental investigation of the dynamic behavior of these microscopic systems is limited by their scale, the identification of macroscale structures capable of describing their response is of particular interest as those offer opportunities to visualize the underlying molecular mechanisms.

In this work, we focus on a chain of pin-joined elastic beams subjected to homogeneous static precompression and use a combination of experiments, numerical simulations and theoretical analyses to investigate the propagation of nonlinear pulses. We find that, while large amplitude rarefaction waves propagate in the chain with constant velocity while conserving their spatial shape, the excited compression waves are strongly dispersive. Further theoretical analysis via a continuum model reveals that the system supports solitary solution only for rarefaction pulses—a behavior that roots from the softening behavior of the beams upon compression. Remarkably, we also find that the model describing our system closely resembles those introduced to characterize the dynamics of several molecular chains and macromolecular crystals [22, 24–26]. As such, since the propagation of pulses in our system can be easily visualized, we envision our model to provide opportunities to elucidate how nonlinear vibrations and pre-deformation affect the macroscopic properties of polymers and other macro-molecular chains.

Our sample consists of a long chain of N = 40 elastic beams free to rotate at their ends, but constrained to move only in longitudinal direction (see figure 1(a)). All beams are made of polyester plastic sheets (Artus Corporation, NJ) and have thickness t = 0.75 mm, length l_b = 50 mm, width b = 60 mm and Young's modulus E = 4.3 GPa. Four pairs of screws and nuts are added along their middle line to increase their mass, therefore reducing the speed of the propagating pulses and facilitating their tracking (see figures 1(a), (b)). Moreover, both ends of the beams are connected via plastic clamps (McMaster-Carr part number 8876T11) to LEGO axles (LEGO part number 3708) covered by sleeves (LEGO part number 62462 and 6590) that serve as hinges. The resulting unit cell has length a = 75 mm and the beam offset by a distance e = 5.2 mm from the the line of action of the applied axial load (see figure 1(b)). Finally, to constrain the chain to move only in longitudinal direction, all LEGO hinges are confined to slide in a metallic rail, lubricated to minimize friction.

Zoom In Zoom Out Reset image size

All our experiments consist of two steps and are conducted on the chain with the left-end connected to a heavy and rigid body and the right-one fixed (see figure 2(c)). First, we compress the structure by slowly moving the heavy rigid body at the left end towards the right by a distance Δx. Such applied displacement induces a longitudinal pre-strain

$$\begin{eqnarray}&&{\varepsilon }_{{\rm{st}}}=-\displaystyle \frac{{\rm{\Delta }}x}{N\,a},\end{eqnarray} \tag{ 1 }$$

in all beams (with only small variations of ∼5% between the units along the chain) and bends them (see figure 2(a)). In figure 2(b) we show the experimentally measured force required to apply a pre-strain ε_st to our beams. We find that the response of the beams is continuous (buckling snap-through is absent due to the eccentricity e) and that their stiffness monotonically decreases as the precompression increases. As such, our pre-deformed elastic units exhibit a strain-softening behavior if further compressed and a strain-hardening response under stretching.

Zoom In Zoom Out Reset image size

Second, we initiate an elastic pulse at the left end of the precompressed chain. We strike the heavy rigid body with a hammer toward the right to excite compression pulses, whereas we simply release the rigid body to generate rarefaction pulses (see figure 2(c)). We then record the elastic wave propagation through the first 30 units of the chain with a digital camera (SONY RX100) at 480 fps (see movie S1 is available online at stacks.iop.org/NJP/21/073008/mmedia) and use digital image correlation [31, 32] to monitor the displacement q_i (see figure 2(c)) of the ith hinge induced by the pulse. Finally, we calculate the longitudinal strain induced by the propagating pulse in the ith beam as

$$\begin{eqnarray}&&{\varepsilon }_{i}=\displaystyle \frac{{q}_{i+1}-{q}_{i}}{{a}_{{\rm{st}}}},\end{eqnarray} \tag{ 2 }$$

where a_st = (1 + _st)a is the unit cell length after precompression.

In figures 2(d)–(g) we focus on two different experiments conducted on a chain pre-compressed by applying ε_st = −0.1 and report both the displacement signal prescribed by the impactor to the first joint, q₁(t), and the spatio-temporal evolution of the strain ε_i. We find that, when we initiate a rarefaction pulse (see figure 2(d)), the wave propagates in a solitary fashion, maintaining both its shape (with an amplitude $\max ({\varepsilon }_{i})\sim 0.11$) and velocity (∼10.4 m s⁻¹—see figure 2(e)). By contrast, when the impact moves the first beam rightwards (see figure 2(f)), the excited compression pulse disperses as it travels through the structure (see figure 2(g)). As such, in full agreement with previous studies on mechanical chains exhibiting strain-softening [16, 33], our experimental results suggest that large amplitude rarefaction waves are stable, whereas compression pulses disperse. Finally, we want to point out that, while the experimental results reported in figures 2(d)–(g) are for a chain with a pre-strain ε_st = −0.1, qualitatively similar behaviors are observed in our tests for ε_st = −0.2 (see figures 2(h), (i)). However, using our setup we could not test propagation of pulses in chains subjected to larger compressive pre-strains, as these result in failure of the beams.

To test the validity of our experimental observations, we establish a discrete model in which each elastic beam is modeled as two rigid rods with length $l=\sqrt{{e}^{2}+{a}^{2}/4}=38\,\mathrm{mm}$ initially rotated by a small angle ${\theta }_{0}=\arctan (2e/a)=0.14$ to account for the eccentricity e (see figure 3(a)). Each pair of rods is connected at the center by a rotational spring with stiffness k_θ = EI/l_b = 0.22 Nm (I = bt³/12 being the area moment of inertia of the beams) that captures the bending stiffness of the elastic beams, while their other two ends are free to slide and rotate (see figure 3(a)). Moreover, we assume that the mass m of the beams is concentrated at the ends of the rigid rods. Specifically, we place a mass m⁽¹⁾ = αm at the ends connected by the torsional spring (with α ∈ [0, 1]) and two masses m⁽²⁾ = (1 − α) m/2 at the other ends (see figure 3(a)). Since in our structure the mass of the beams, hinges and screws/nuts is 3 g, 7 g and 8 g, we find that for the considered system m = 18 g and α = 0.56. It should also be noted that our discrete model is similar to the one proposed to describe solitons propagation in polyethylene macromolecules [24–26], providing an interesting analogy between mechanical beams and molecular units.

Zoom In Zoom Out Reset image size

Identically to our experiments, all our simulations consists of two steps. In the first step we statically compress the structure by reducing the unit cell length to a(1 + ε_st) (see figure 3(a)). Such reduction in unit cell length increases the angle between all rods and the horizontal line from θ₀ to

$$\begin{eqnarray}&&{\rm{\Theta }}={\theta }_{0}+{\theta }_{{\rm{st}}},\end{eqnarray} \tag{ 3 }$$

where θ_st is the rotation induced by the precompression, which can be expressed as a function of ε_st using

$$\begin{eqnarray}&&\cos {\rm{\Theta }}=\displaystyle \frac{a}{2l}(1+{\varepsilon }_{{\rm{st}}}).\end{eqnarray} \tag{ 4 }$$

Note that the precompression of the beams requires application of a force F_st that satisfies

$$\begin{eqnarray}&&{F}_{{\rm{st}}}\,l\sin {\rm{\Theta }}-{k}_{\theta }{\theta }_{{\rm{st}}}=0.\end{eqnarray} \tag{ 5 }$$

In figure 3(b) we compare the F_st–ε_st relation predicted by equations (3)–(5) with experimental measurements and find that the static response of the beams is well captured by the discrete model.

In the second step, we then simulate the propagation of nonlinear waves in the pre-deformed chain (see figure 3(c)). To this end, we define the Lagrangian of the system

$$\begin{eqnarray}&&L=\displaystyle \sum _{i=1}^{N}\ \displaystyle \frac{1}{2}{m}^{(1)}\left({\dot{u}}_{i}^{2}+{\dot{{v}_{i}}}^{2}\right)+\displaystyle \frac{1}{2}{m}^{(2)}{\dot{q}}_{i}^{2}-\displaystyle \frac{1}{2}{k}_{\theta }{\left(2{\theta }_{{\rm{st}}}+2{\theta }_{i}\right)}^{2},\end{eqnarray} \tag{ 6 }$$

where u_i and v_i are the longitudinal and transverse displacements of the ith mass m⁽¹⁾ excited by the pulse and θ_i denotes the wave-induced rotation of the ith pair of rods (see figure 3(c)). Since in our discrete model each unit has only one degree of freedom, u_i, v_i and θ_i can be expressed as a function of the longitudinal displacement of the hinges, q_i, as

$$\begin{eqnarray}&&{u}_{i}=\displaystyle \frac{1}{2}({q}_{i+1}+{q}_{i}),\end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray}&&{v}_{i}=\sqrt{{l}^{2}-{\left(\displaystyle \frac{{q}_{i+1}-{q}_{i}+2l\cos {\rm{\Theta }}}{2}\right)}^{2}}-l\sin {\rm{\Theta }},\end{eqnarray} \tag{ 7b }$$

$$\begin{eqnarray}&&{\theta }_{i}=\arccos \left(\displaystyle \frac{{q}_{i+1}-{q}_{i}}{2l}+\cos {\rm{\Theta }}\right)-{\rm{\Theta }}.\end{eqnarray} \tag{ 7c }$$

By using equations (7), (6) can be written only in terms of q_i and ${\dot{q}}_{i}$, and the discrete motion equations of the system can then be obtained via the Euler–Lagrange equations as

$$\begin{eqnarray}&&\displaystyle \frac{{\rm{d}}}{{\rm{d}}t}\left[\displaystyle \frac{\partial L({q}_{i},{\dot{q}}_{i})}{\partial {\dot{q}}_{i}}\right]-\displaystyle \frac{\partial L({q}_{i},{\dot{q}}_{i})}{\partial {q}_{i}}=0,\ \ i=1,\,..,N.\end{eqnarray} \tag{ 8 }$$

For a chain comprising N beams, equation (8) result in a system of N coupled differential equations, which we numerically solve using the 4th order Runge–Kutta method (via the Matlab function ode45—see supporting information for the Matlab code). To test the relevance of our discrete model, in figures 3(d) and (e) we focus on the two tests presented in figure 2 and compare the evolution of ε_i along the chain as extracted from experiments (square markers) and simulations (triangular markers) at three different times. Moreover, in figures 3(f) and (g) we show the numerically obtained phase-space plots for the two pulses. Note that in our numerical analysis we consider a chain comprising N = 40 units, assign ε_st = −0.1 to all beams (as in our experiments), apply the experimentally extracted displacement signal q₁(t) (see figures 2(d) and (f)) to the first beam and implement fixed boundary conditions at the right end. Remarkably, the numerical results nicely capture the behavior observed in our experiments. As shown in figure 3(d), the numerical analyses indicate that the excited rarefaction pulse propagates without apparent distortion with a velocity of ∼9.8 m s⁻¹. This prediction is very close to the experimentally measured value of ∼10.4 m s⁻¹—with the discrepancy mostly arising because our simple model underestimates the stiffness of the beams upon loading (see figure 3(a)). Moreover, the phase-space plot for rarefaction shows a clear homoclinic orbit (see figure 3(f)), which is consistent with a solitary wave solution. By contrast, the compressive pulse is not stable (see figure 3(e))—a feature that is also captured by the chaotic trajectories emerging in the phase-space plot (see figure 3(g)).

Having verified the validity of our discrete model, we then simplify equations (8) to derive an analytical solution. To this end, we first introduce a continuous function q that interpolates the discrete variables q_i as

$$\begin{eqnarray}&&q(x={x}_{i})={q}_{i},\end{eqnarray} \tag{ 9 }$$

where x_i = (i − 1)a_st denotes the position of the left end of the ith beam in the chain after static precompression. Then, we assume that: (I) ${\varepsilon }_{i}\ll 1$, so that equations (7b) and  (7c) can be approximated as

$$\begin{eqnarray}\begin{array}{rcl}{v}_{i} & \approx & -\displaystyle \frac{{q}_{i+1}-{q}_{i}}{2\tan {\rm{\Theta }}}-\displaystyle \frac{{({q}_{i+1}-{q}_{i})}^{2}}{8l{\sin }^{3}{\rm{\Theta }}},\\ {\theta }_{i} & \approx & -\displaystyle \frac{{q}_{i+1}-{q}_{i}}{2l\sin {\rm{\Theta }}}-\displaystyle \frac{\cos {\rm{\Theta }}\,{({q}_{i+1}-{q}_{i})}^{2}}{8{l}^{2}{\sin }^{3}{\rm{\Theta }}},\end{array}\end{eqnarray} \tag{ 10 }$$

and (II) the width of the propagating pulse is much larger than the unit cell length a, so that the displacement of the (i − 1)th and (i + 1)th joints can be expressed using Taylor expansion as

$$\begin{eqnarray}\begin{array}{rcl}{q}_{i+1} & \approx & {\left.\left(q+{a}_{{\rm{st}}}{q}_{x}+\displaystyle \frac{{a}_{{\rm{st}}}^{2}}{2}{q}_{{xx}}+\displaystyle \frac{{a}_{{\rm{st}}}^{3}}{6}{q}_{{xxx}}+\displaystyle \frac{{a}_{{\rm{st}}}^{4}}{24}{q}_{{xxxx}}\right)\right|}_{x={x}_{i}}\\ {q}_{i-1} & \approx & {\left.\left(q-{a}_{{\rm{st}}}{q}_{x}+\displaystyle \frac{{a}_{{\rm{st}}}^{2}}{2}{q}_{{xx}}-\displaystyle \frac{{a}_{{\rm{st}}}^{3}}{6}{q}_{{xxx}}+\displaystyle \frac{{a}_{{\rm{st}}}^{4}}{24}{q}_{{xxxx}}\right)\right|}_{x={x}_{i}},\end{array}\end{eqnarray} \tag{ 11 }$$

where q_x = ∂q/∂x. Substitution of equations (10), (11) into equation (8) yields the continuum governing equation

$$\begin{eqnarray}\begin{array}{rcl}{q}_{{tt}} & = & {c}_{0}^{2}(1-{\theta }_{{\rm{st}}}\cot {\rm{\Theta }})\left({q}_{{xx}}+\displaystyle \frac{{l}^{2}}{12}\ {q}_{{xxxx}}\right)+\alpha {l}^{2}({\cot }^{2}{\rm{\Theta }}-1){q}_{{ttxx}}/4\\ & & +{c}_{0}^{2}\left(3\cot {\rm{\Theta }}-{\theta }_{{\rm{st}}}-3{\theta }_{{\rm{st}}}{\cot }^{2}{\rm{\Theta }}\right)\cot {\rm{\Theta }}\,{q}_{x}{q}_{{xx}},\end{array}\end{eqnarray} \tag{ 12 }$$

where ${c}_{0}=2{a}_{{\rm{st}}}\cot {\rm{\Theta }}\sqrt{{k}_{\theta }/m}$ denotes the characteristic velocity of the system. Finally, we take the derivative of equation (12) with respect to x and then introduce the continuum strain distribution ε = ∂q/∂x to obtain

$$\begin{eqnarray}\begin{array}{rcl}{\varepsilon }_{{tt}} & = & {c}_{0}^{2}(1-{\theta }_{{\rm{st}}}\cot {\rm{\Theta }})\left({\varepsilon }_{{xx}}+\displaystyle \frac{{l}^{2}}{12}\ {\varepsilon }_{{xxxx}}\right)+\alpha {l}^{2}({\cot }^{2}{\rm{\Theta }}-1){\varepsilon }_{{ttxx}}/4\\ & & +\displaystyle \frac{{c}_{0}^{2}}{2}\left(3\cot {\rm{\Theta }}-{\theta }_{{\rm{st}}}-3{\theta }_{{\rm{st}}}{\cot }^{2}{\rm{\Theta }}\right)\cot {\rm{\Theta }}\,{\left({\varepsilon }^{2}\right)}_{{xx}},\end{array}\end{eqnarray} \tag{ 13 }$$

which has the form of a double-dispersion Boussinesq equation [34] with the double-dispersion term ε_ttxx introduced because of the intertia coupling between m⁽¹⁾ and m⁽²⁾ (if m⁽¹⁾ = 0, then α = 0 and the coefficient in front of ε_ttxx vanishes). Note that the double-dispersion Boussinesq equation was first derived by Boussinesq in the 19th century to account for both the horizontal and vertical flow velocity in the description of nonlinear water waves [34] and has subsequently been used to describe nonlinear waves in a variety of systems, including microstructured solids [35, 36] and polyethylene [24].

Remarkably, it has been shown that the Boussinesq equation admits the solitary wave solution [37, 38]

$$\begin{eqnarray}&&\varepsilon =A\,{{\rm{{\rm{sech}} }}}^{2}\left(\displaystyle \frac{x-ct}{W}\right),\end{eqnarray} \tag{ 14 }$$

where c is the velocity of the pulse and A and W are its amplitude and characteristic width, which are given by

$$\begin{eqnarray}\begin{array}{rcl}A & = & \displaystyle \frac{3\left({c}^{2}-{c}_{0}^{2}+{c}_{0}^{2}\ {\theta }_{{\rm{st}}}\cot {\rm{\Theta }}\right)}{{c}_{0}^{2}(3\cot {\rm{\Theta }}-{\theta }_{{\rm{st}}}-3{\theta }_{{\rm{st}}}{\cot }^{2}{\rm{\Theta }})\cot {\rm{\Theta }}},\\ W & = & l\sqrt{\displaystyle \frac{\tfrac{{c}_{0}^{2}}{3}\left(1-{\theta }_{{\rm{st}}}\cot {\rm{\Theta }}\right)+\alpha \left({\cot }^{2}{\rm{\Theta }}-1\right){c}^{2}}{{c}^{2}-{c}_{0}^{2}+{c}_{0}^{2}\ {\theta }_{{\rm{st}}}\cot {\rm{\Theta }}}}.\end{array}\end{eqnarray} \tag{ 15 }$$

Finally, the displacement distribution q(x) can be obtained as

$$\begin{eqnarray}&&q(x)=\int \varepsilon (x){\rm{d}}x={AW}\tanh \left(\displaystyle \frac{x-ct}{W}\right)+C\end{eqnarray} \tag{ 16 }$$

where C is the integration constant that is found as C = −AW by imposing $q(x\to \infty )=0$.

Using equation (15), we find that for A = 0.11 (i.e. for a pulse with the same amplitude as the one excited in the experiment of figures 2(d) and (e)), a solitary wave with c = 9.31 m s⁻¹ and W = 91.6 mm propagates through the system, in good agreement with both our experimental and numerical observations (see figure 3(d)). By contrast, for a compression pulse with A = −0.11 (i.e. for a pulse with the same amplitude as the one excited in the experiment of figures 2(f) and (g)), equation (15) yields an imaginary W, confirming that our solitary wave solution is not valid anymore.

Having verified the ability of the analytical solutions to capture our experimental observations, we now use it to investigate how the pre-strain ε_st affects the propagation of the pulses. In figures 4(a) and (b) we report the evolution of W and c predicted by equation (15) as a function of A and ε_st. Remarkably, we find that, irrespective of ε_st, our structure supports solitons only if A > 0. For A < 0 the width W predicted by equation (15) is an imaginary number and our solitary solution is not anymore valid. We also find that the width W of the supported rarefaction solitons is mostly affected by the amplitude A, whereas the pulse velocity c monotonically decreases with the increase in precompression (−ε_st)—a dependency that is consistent with the strain-softening behavior of the unit cells under compression. Finally, we note that the predictions of our continuum model are in excellent agreement with those obtained by directly integrating the system of ordinary differential equations given by equation (8). As an example, in figures 4(c)–(h) we compare the results given by our continuum and discrete models for A = ±0.11 and three different levels of applied pre-strain (i.e. ε_st = −0.1, −0.2 and −0.3). The numerical results confirm that only for A > 0 the pulses are stable (and characterized by width and speed very close to those predicted by equation (15)), whereas for A < 0 they are strongly dispersive.

Zoom In Zoom Out Reset image size

Lastly, we use our analytical solution to investigate the influence of the structural parameters (α, m, k_θ and θ₀) on the characteristics of the supported nonlinear waves. We find that, while m and k_θ (which can be altered by changing the dimensions and material properties of the beams) only affect the characteristic velocity of the system ${c}_{0}=2{a}_{{\rm{st}}}\cot {\rm{\Theta }}\sqrt{{k}_{\theta }/m}$, α (which can be altered by adding or removing screws and nuts along the center line of the beams) and θ₀ (which can be tuned by changing the offset distance e) alter both c and W (see figure 5). Furthermore, an exhaustive search in all combinations of structural parameters predicts that nonlinear pulses are stable if and only if A > 0 (see figure 5)—indicating that this system can support rarefaction solitons but not compression ones.

Zoom In Zoom Out Reset image size

To summarize, we use a combination of experiments, numerical analyses and theory to investigate the propagation of solitary waves in a 1D chain of precompressed elastic beams. First, we have conducted experiments on a centimeter-scale system to characterize the propagation of large amplitude compression and rarefaction pulses and found that, while the compression pulses disperse, the rarefaction ones retain their shape and propagate with constant velocity. Second, we have derived a continuum model that captures the experimental observations and confirms that the system supports rarefaction solitons only. It is important to emphasize that such behavior is consistent with the static behavior of the precompressed beams (see figure 3(b)), since it has been shown that a chain of units exhibiting softening behavior only supports the propagation of rarefaction solitary waves [33].

This study represents the first step towards the investigation of large amplitude waves in beam lattices in 2D and 3D. While periodic lattices have recently attracted considerable interest because of their ability to tailor the propagation of linear elastic waves through directional transmissions and band gaps (frequency ranges of strong wave attenuation) [39–43], comparatively little is known about their nonlinear behaviors under high-amplitude impacts [20]. The results presented in this paper provide useful guidelines for future explorations of the propagation of nonlinear waves in lattice materials. Finally, since the derived model shares strong similarities with those established to describe the dynamics of polymer and macromolecular crystals [22–27], we believe that our experimental platform (which gives direct access to all parameters and variables of the system) could provide insights into a range of nonlinear wave effects relevant at the molecular scale. In particular, the extension of the present system to a bistable one in which the beams could bend up and down (in our current setup the bending direction is limited by the surface on which the system is placed) could elucidate the dynamics of topological solitons in polyacetylene [23].

KB acknowledges support from the National Science Foundation under Grant No. DMR-1420570 and EFMA-1741685 and from the Army Research Office under Grant No. W911NF-17-1-0147.