Computational model of twisted elastic ribbons

We develop an irregular lattice mass-spring model to simulate and study the deformation modes of a thin elastic ribbon as a function of applied end-to-end twist and tension. Our simulations reproduce all reported experimentally observed modes, including transitions from helicoids to longitudinal wrinkles, creased helicoids and loops with self-contact, and transverse wrinkles to accordion self-folds. Our simulations also show that the twist angles at which the primary longitudinal and transverse wrinkles appear are well described by various analyses of the Föppl–von Kármán equations, but the characteristic wavelength of the longitudinal wrinkles has a more complex relationship to applied tension than previously estimated. The clamped edges are shown to suppress longitudinal wrinkling over a distance set by the applied tension and the ribbon width, but otherwise have no apparent effect on measured wavelength. Further, by analyzing the stress profile, we find that longitudinal wrinkling does not completely alleviate compression, but caps the magnitude of the compression. Nonetheless, the width over which wrinkles form is observed to be wider than the near-threshold analysis predictions: the width is more consistent with the predictions of far-from-threshold analysis. However, the end-to-end contraction of the ribbon as a function of twist is found to more closely follow the corresponding near-threshold prediction as tension in the ribbon is increased, in contrast to the expectations of far-from-threshold analysis. These results point to the need for further theoretical analysis of this rich thin elastic system, guided by our physically robust and intuitive simulation model.

Figure 1

Deformation modes of twisted thin sheets, replicated by our simulation whose details are summarized in Fig. 2. Ribbon (a) is labeled with its dimensions and quantities relevant to the twisting-under-tension procedure. $T$ is the scaled longitudinal tension applied, which is total applied force $F$ normalized by the width, Young's modulus, and thickness of the ribbon, and $\eta$ is a scaled twist angle: the end-to-end angle $\theta$ normalized by the ribbon's aspect ratio. These quantities are defined in Table 1. (a) The helicoid phase initially present before any buckling transitions occur. (b) Longitudinal buckling occurs at an angle ${\eta }_{\text{lon}}$ and has a fixed wavelength ${\lambda }_{\text{lon}}$ set by the tension $T$ and thickness $h$. (c) Creased helicoids develop from the longitudinally buckled ribbons as the buckle ridges “turn” to form triangular facets. (d) At tensions below the crossover tension $T^{*}$, ribbons will snap through to form a loop at an angle ${\eta }_{\text{tran}}$. If twisted far enough, ribbons will develop self-contact at an angle ${\eta }_{\text{sc}}$. (e) Transverse buckling occurs at tensions greater than $T^{*}$, transitioning at angle ${\eta }_{\text{tran}}$. (f) The wavelength of transverse buckling ${\lambda }_{\text{tran}}$ is set by the aspect ratio of the sheet; sheets with smaller aspect ratios can display several wavelengths of transverse buckling. (g) At high twists, low aspect ratio sheets display “accordion folding” and approach the yarning transition [28]. Snapshots (a)–(e) fall within the phase diagram presented in Fig. 3.

Figure 2

Nodes are connected by in-plane springs with stiffness $k_s$ and dashpot damping $b_{\text{int}}$. The $k_s$ of each in-plane spring is set by the local geometry of the mesh and the target Young's modulus for the sheet, $Y$. Out-of-plane stiffness $k_b$ is similarly dictated by the local geometry and the target bending rigidity, $B$. Both sets of springs ensure quadratic energy penalization for stretching and misalignment of the vectors normal to the triangular mesh facets, Eqs. (4) and (5).

Figure 3

Slice of the thin ribbon deformation phase space at fixed thickness $h=127\mu \mathrm{m}$. Plotted are longitudinal buckling transitions (blue circles, long-dashed border), transverse buckling transitions (green circles, short-dashed border), and points of self-contact (red circles, solid border) as a function of scaled ribbon tension $T$ and scaled twist angle $\eta$. At low tensions the transverse buckling transition often leads immediately to self-contact, whereas at higher tensions the transverse instability occurs before self-contact develops. The blue long-dashed line corresponds to the theoretical scaling of the longitudinal buckling angle ${\eta }_{\text{lon}}$ [Eq. (14)] with a single fitting parameter for the coefficient of the finite thickness correction term. The green short-dashed line is the theoretical scaling of the transverse buckling angle ${\eta }_{\text{tran}}$ [Eq. (15)], again with a single fitting parameter for the coefficient. $T^{*}$, the tension at which the primary instability switches from longitudinal to transverse buckling, is indicated by the vertical black dashed line.

Figure 4

Detection of longitudinal wrinkles. (a) The out-of-plane displacement $\widehat{z}$ along the midline of a sample ribbon $127\mu \mathrm{m}$ thick with scaled tension $T=9.00\times {10}^{-4}$, which undergoes longitudinal buckling. Note that due to the clamped boundary conditions the wrinkles are suppressed a distance $L_{\text{sup}}\sim O\left(W\right)$. (b) The amplitude spectrum $a\left(f\right)$ of the middle third of the wrinkle profile in (a), computed as the absolute value of the discrete Fourier transform of $\widehat{z}$. The frequency of the highest peak is taken to be the reciprocal of the wrinkle wavelength $\lambda$. To aid in identifying the precise onset of wrinkling, the amplitude spectrum is truncated to retain only frequencies in the range $(1/\lambda )\pm \mathrm{\Delta }$, and the wrinkle profile is reconstructed from the inverse Fourier transform of the truncated spectrum. (c) A waterfall plot of reconstructed wrinkle profiles as a function of scaled twist angle $\eta$ along the $y$ axis, for which frequencies $(1/\lambda )\pm 0.15\mathrm{c}{\mathrm{m}}^{-1}$ have been retained, with $\lambda \approx 1.69\mathrm{c}\mathrm{m}$. (d) The wrinkle amplitude ${\left[{\langle H^2(r=0)\rangle }_y\right]}^{1/2}$ as a function of $\eta$, where the average ${\langle ·\rangle }_y$ is computed over the middle third of the wrinkle profiles in (c). The three identified scaled angles ${\eta }_{\text{lon}}$, ${\eta }_1$, and ${\eta }_2$ correspond to selected wrinkle profiles in (c) shown in black. The scaled angle ${\eta }_{\text{lon}}=0.203$ marks the detected onset of longitudinal wrinkling, identified as the “knee point” [34] of the amplitude curve. Inset: Amplitude as a function of confinement parameter offset by the wrinkle onset, $\alpha -{\alpha }_{\text{lon}}$, with scaling $\beta =0.56$ [Eq. (17)].

Figure 5

Detection of transverse buckling and self-contact. (a) The shortest distance between any two points on opposing long edges of a ribbon as a function of scaled twist angle, shown for ribbons at three different scaled tensions. At $\eta =0$, the edge distance is equal to the ribbon width $W$, while after self-contact, the edge distance equals the effective thickness $h$. Solid markers indicate the onset of longitudinal buckling, transverse buckling, and self-contact, mapping out a vertical trajectory through the phase space plot of Fig. 3 along constant $T$. (b) Snapshots of the ribbon in (a) (middle) at $T=1.98\times {10}^{-3}$ in each of the four deformation regimes.

Figure 6

Variation of scaled longitudinal wavelength ${\lambda }_{\text{lon}}/\sqrt{hW}$ with scaled tension $T$. (a) The wavelength ${\lambda }_{\text{lon}}$ of the wrinkles was estimated in experiments to scale as $T^{-1/4}$, and good agreement is obtained by fitting a curve through the data at all but the two smallest tensions ($T<4\times {10}^{-4}$), as shown by the dashed line [9]. Prior experimental results also reveal deviations from the theoretical scaling at low tensions [9]. (b) Individual fits of scaled wavelength for each of the three thickness to width ratios $h/W$ (solid lines) shown on a log-log scale, with the theoretical scaling reproduced from (a) (dashed line). Closer examination of the data from (a) reveals that ${\lambda }_{\text{lon}}$ has a lingering dependence on the thickness $h$, and could have a more complex relationship to $T$ than previously estimated.

Figure 7

Wrinkle suppression in a ribbon with $h=127\mu \mathrm{m}$ due to clamped boundaries. (a) Profiles of all longitudinally buckled samples, highlighting the wrinkle suppression zone near the ribbon edges in blue. The suppression distance measured from either clamp is $L_{\text{sup}}$. (b) Wrinkle suppression length, $L_{\text{sup}}$ as a function of applied tension $T$. The suppression near a clamped edge approaches the value $2W=5.08\mathrm{c}\mathrm{m}$ as tension increases toward $T^{*}$. $L_{\text{sup}}$ likely depends on $h$ as well as $W$.

Figure 8

(a) Stress profile of a ribbon with $T=1.08\times {10}^{-3}$, $h=127\mu \mathrm{m}$. The color bar (and color of each stress slice) indicates the progression of the ribbon's twist, with blue corresponding to $\eta =0$ and red $\eta =0.8$, which is the fully twisted state for this sample. The solid black line in the color bar marks ${\eta }_0$, when stress becomes compressive, and the gray dashed line marks ${\eta }_{\text{lon}}$, the onset of longitudinal buckling. As the ribbon twists, its longitudinal stress profile is initially parabolic, crossing over to compressive stress in the center of the ribbon at ${\eta }_0$ ($=0.161$). After the longitudinal buckling transition at ${\eta }_{\text{lon}}$ ($=0.216$, dashed blue-purple line) the stress profile bottoms out (minimum stress marked by the horizontal gray, dashed line) and widens. The buckled ribbon thus continues to support compressive stress but functionally caps the compressive load. (b) At varying thicknesses ($h=63.5\mu \mathrm{m},127\mu \mathrm{m},254\mu \mathrm{m},508\mu \mathrm{m}$) and degrees of confinement ($\alpha =80,120,160,200$), we see that the residual stress supported by the longitudinally buckled ribbon is linear in $h/\left(W\sqrt{T}\right)$, with a slope of $-5.14\pm 0.05$ given by the dashed line. Regardless of $\alpha$, as the thickness tends to zero we expect the stress in the ribbon to also vanish.

Figure 9

Development of the wrinkles that appear in a ribbon with a tension $T=1.08\times {10}^{-3}$, $h=127\mu \mathrm{m}$. Snapshots were taken at regular intervals, starting at the onset of buckling. The first frame is at $\eta ={\eta }_{\text{lon}}$, and the frames below show the progression of the wrinkles as the ribbon goes through the longitudinally buckled phase and begins to cross into the creased helicoid phase. (a) Contour maps of the mean curvature in the ribbon. Red (blue) corresponds to above(below)-the-plane curvature. Black (gray) dashed lines mark the near threshold (far-from-threshold) predictions for the width of the wrinkles. At the onset of wrinkling their width is consistent with the NT prediction, but they soon expand to meet the FT prediction. In the final frame we also see the wrinkle ridges begin to “turn” and start to form the triangular facets characteristic of the creased helicoid phase. (b) Contour maps of the longitudinal stress component. Red (blue) corresponds to tensile (compressive) in-plane stress. As the twist angle increases, the magnitude of the compressive stress remains approximately the same, though the width of the compressive area broadens. The transition from longitudinal buckling to the angled creased helicoid is also evident in these stress maps. (c) Stress profile and resultant amplitude $[A={\left({\langle H^2\left(r\right)\rangle }_y\right)}^{1/2}]$ at various values of $\eta$.

Figure 10

Maximum amplitude of longitudinal wrinkles scales like $A(r=0)\propto {(\alpha -{\alpha }_{\text{lon}})}^{\beta }$. Values of $\beta$ are extracted from the amplitude curve in Fig. 4 and are plotted here for ribbons of varying thickness and applied tension. We find $\beta =0.59\pm 0.05$ with the average value given by the solid black line. The horizontal dashed black line is the value Chopin and Kudrolli [12] experimentally estimated to be $\beta =2/3$.

Figure 11

Tests of the NT (red, upper lines) and FT (blue, lower lines) model predictions at $h=127\mu \mathrm{m}$, various ribbon tensions. (a) At the onset of wrinkling (${\eta }_{\text{lon}}$, dashed black line) the wrinkled zone emerges at a nonzero width. It widens as the twist progresses, eventually asymptoting to $r_{\text{wr}}=0.5$ as $\alpha \to \infty$. This behavior is similar across various applied tensions, though at moderate tensions (e.g., $T=1.08\times {10}^{-3}$) $r_{\text{wr}}$ adheres more closely to the theoretical FT prediction than it does at lower tensions. (b) Another way to test the NT and FT predictions is by measuring the contraction $\chi$ of the ribbon. Here $\chi$ is plotted against the inverse confinement parameter (such that small values of $1/\alpha$ correspond to the greatest angles $\eta$). Across varying tensions the ribbon tends to contract as a pure helicoid, following the NT prediction, until the onset of wrinkling (${\eta }_{\text{lon}}$, dashed black line). After wrinkling, the contraction proceeds rapidly, at times even faster than the predicted FT slope. This difference in scaling could be because of the transition to the creased helicoid phase, which is predicted to have a different contraction scaling than the FT model [7, 13, 14]. The gray dashed line corresponds to the transverse buckling point, after which the ribbon springs back to a lesser contraction.

Figure 12

Uniaxial stretching errors in the regular (left, triangular markers) and random (right, X markers) meshes as a function of mesh resolution. For both meshes the mVG model (blue, lower point clouds) error has a smaller magnitude than the SN model (red, upper point clouds).

Figure 13

Percent errors in the measured Poisson ratio compared to the expected value of $\nu =1/3$ in the regular (left, triangular markers) and random (right, cross markers) meshes. The mVG model (blue, lower point clouds) more reliably produces the expected Poisson ratio than the SN model (red, upper point clouds).

Figure 14

Pure shearing errors in the regular (left, triangular markers) and random (right, cross markers) meshes as a function of mesh resolution. For both meshes the SN model (red, lower point clouds) has smaller magnitude of errors and converges at a higher order than the mVG model (blue, upper point clouds).

Figure 15

Uniformly loaded, simply supported plate bending errors in the regular (left, triangular markers) and random (right, X markers) meshes as a function of mesh resolution. For the regular mesh the mG model (blue, lower point clouds) has a smaller magnitude of error and greater order of convergence than the SN model (red, upper point clouds). The random meshes converge at the same order, but in general the magnitude of the error is smaller when using the mG model.