Tensional twist-folding is a method to transform flat sheets into layered structures and yarns with ordered internal architectures by remote boundary manipulation. Twisting sheets under tension has been used since antiquity in making catgut bow strings, surgical sutures, musical chord instruments, sports rackets, sausage and candy wrappers, fabric filters and wearbles such as turbans and crushed dupattas, and in upcycling of plastic (Fig. 1 A to C). Scrolled yarns with nested structures optimized for energy harnessing, batteries, and embedding materials (1)(2)(3)(4) are difficult to achieve by compression-induced transformations of elastic sheets (5)(6)(7), and traditional fiber spinning methods (8). When appropriate materials are used, the transformations can be reversible, and twist folding and scrolling can be used to reconfigure and repurpose flat sheets, as exemplified by the multifunctional Rajastani turban.
The interplay between topology and large shape transformations have been studied in terms of inextensible twisted rods and ribbons to understand the conformation of DNA and proteins (9)(10)(11), and has contributed to the development of a now well-established theoretical framework (12,13). However, shape transformation of sheets which significantly stretch upon twist have remained undocumented despite their ubiquity in a wide range of applications.
We report the spontaneous formation of twisted multi-layered yarns with ordered internal architectures enabled by x-ray 3D scanning. These structures obtained under extreme deformation From [25] PDMS, t/W = 0.02 and self-contact are distinct from those observed in rods and ribbons at moderate twist reported previously (14,15), and are not known to occur by a purely compression-driven transformations of elastic sheets as in crumpling, folding, and capillary wrapping (16)(17)(18)(19). Modeling such large shape transformations and configurations is extremely challenging. Elastic plate models such as the Föppl-von Kármán (FvK) equation, and its more recent co-variant extension (20), have solved the initial growth of wrinkling above onset of primary instability (20)(21)(22)(23)(24)(25)(26), but fail to anticipate, let alone explain, the proposed transformation of a flat sheet into scrolled yarns with functional guests (2). As in other paradigm pattern formation systems such as buoyancydriven Rayleigh-Bénard convection which displays intermittent spatio-temporal chaos (27), it is not a priori obvious what imprints of the primary instabilities persist as a sheet is twisted far beyond the perturbative regime where previous studies were focused. While origami and inextensible sheet models are amenable to address large shape transformations (18,(28)(29)(30)(31)(32)(33)(34)(35)(36)(37)(38), their generalization to significantly stretched sheets is unknown.
Going beyond reporting the discovery of ordered transformations, we develop a framework which combines the kinematics of stretched sheets, origami, and fold-induced transverse stiffness to explain our observations. Remarkably, we find that the observed accordion folded sheets have regular polygonal shapes described by Schläfli symbols (39), and show that origami kinematics can capture the main features of the structure. We provide an analytical framework to address the successive transformations experienced by a twisted sheet from the onset of transverse wrinkling via recursive folding and scrolling. Our framework can serve as a guide for fabrication of yarns with precise control of crosssectional architecture. When made with hyperelastic materials, which recover their unstressed states, they can be repeatedly reconfigured and redeployed with our twist-folding method.
Ordered shape-transformation and nonmonotonic torque with twist
Examples of a polydimethylsiloxane (PDMS) sheet with increasing twist are shown in Fig. 1D-G and SM: Movie S1. A system schematic and the Cartesian coordinates system (x, ȳ, z) are shown in Fig. 1D, inset. The system consists of a sheet of length L, width W , and thickness t,
twisted by an angle θ while being held at opposite ends and stretched axially by ∆L. Transverse wrinkles can be observed just above the onset of primary instability (Fig. 1D), which grow in amplitude and collapse into an accordion folded spiral structure with self-contact (Fig. 1E). As θ is increased further, a nestled helical structure forms at the waist (Fig. 1F), before a secondary instability occurs which leads to recursive folding and a scrolled yarn (Fig. 1G). Each of the major shape transformations causes the rate of change of applied torque M to change sign, leading to a sawtooth variation with twist (Fig. 1H). The primary instability and parameter space over which these transformations occur varies with L/W (Fig. 1I). The observed angle at which the primary instability occurs in Fig. 1I is consistent with the clamp-dominated (L ∼ W ) and the ribbon (L W ) regimes, which scales as θ p ∼ (L/W ) ζ , with ζ = 1 and 2, respectively, and wavelength λ p ∼ √ Lt(∆L/L) -1/4 reported previously (20,25). The ordered sheet transformations and sawtooth torque variations reported in Fig. 1E-I as self-contact occurs are documented here for the first time, and are the focus of the following analysis.
We investigate the morphologies of the sheets beyond incipient wrinkling until self-contact, with noninvasive 3D x-ray tomography. In Fig. 2A, we reconstruct the central 80% of a twisted poly-vinyl siloxane (PVS) sheet. We calculate its mean curvature H (see SM 6: Curvature calculation) and note high-curvature regions along wrinkle antinodes which develop above the onset of transverse instability. To recognize the spatial distribution of the curvature as the sheet wraps around itself, we map H to a rectangular domain with axes (s, x), the curvilinear and vertical coordinates, respectively (Fig. 2B, and SM: Figs. S2 andS5A). The wrinkles are observed to be initially aligned with the applied tension when θ = 90 • consistent with linear perturbation analysis (20). With increasing twist, H is increasingly localized along folds with essentially flat regions in between, and the folds rotate away from the tensional axis till they meet near the clamped edges.
We calculate the bending energy density
] is the bending stiffness, ν the Poisson ratio, K the Gaussian curvature, and Y the Young's modulus. Plotting w b across the sheet at mid-distance between the clamps, we observe that sharp peaks grow with θ showing that the bending energy increasingly localizes along the folds (Fig. 2C). Since the number of peaks is unchanged as twist is increased, we postulate that the number of folds n are set by twice the ratio of W and λ p . Accordingly, we plot the observed n versus the predicted n = 2W/λ p in Fig. 2D, and observe good agreement, showing the imprint of the primary instability far into the shape transformation.
Based on these observations, we develop a model to compute the torsional response and stored elastic energy of the sheet as it folds until reaching self-contact. We assume that the sheet folding is predominantly given by elongation along straight lines, as represented by the red and blue lines shown in Fig. 3A connecting the clamped edges and the central crosssection. The elastic energy E el of the sheets is calculated by integration of the strain energy density Y /2 2 where = (L -L)/L is the strain and L the strip length (see also SM 7: Pre-self-contact elastic energy and torque analytics):
The evolution of the strip length with twist may seem simple, but in fact it encodes a complex folding kinematics. Here, we assume that these kinematics can be captured by the compaction parameter χ = D/W , where D is the diameter of a circle enclosing the central sheet crosssection (Fig. 3A and SM: Fig. S6). Then,
where,
which is related to the twist induced strain. Using Eqs. [2] and [3], E el is calculated by numerically integrating Eq. [1], and an analytical expression is obtained as a polynomial of F accurate within 10% error for L > 2W (see Eq. S16 in SM: 7. Pre-self-contact elastic energy and torque analytics).
Recalling that M = dE el /dθ, we can then obtain an analytical expression for the torque:
where M 0 = Y W 2 t. Fig. 3B shows a comparison of the measured M (θ) with that calculated using Eq. [4]. The elasto-geometric fold model is seen to capture the nonmonotonic torsional response quantitatively until the onset of self-contact at θ ≈ 180 • . We show that the nonmonotonicity originates from finite rotation effects which are essentially handled in Eq. [4] by a sine function. Furthermore, the occurrence of torque peak at θ ≈ 105 • is also captured. Because the peak occurs at θ > 90 • , the quantity χ encodes the fold-induced resistance to transverse displacement according to our elasto-geometric fold model. If this bending-induced stiffness was absent, the peak would occur at θ = 90 • (see SM: Fig. S6B). It is noteworthy that M calculated using the perturbative Föppl-von Kármán equations (14,40) predicts an ever increasing torque with twist which is clearly not the case in our data (see SM 7: Pre-self-contact elastic energy and torque analytics, C and D).
We calculate the elastic energy obtained by numerical integration of Eq. [1] as a function of θ, and compare it with the measured values in Fig. 3C after normalizing with the zero twist stretching energy E 0 . We find very good agreement showing that the bending energy can be neglected (Fig. 3D and SM 8: Comparison Stretching and Bending Energy up to Half-turn Twist). It is noteworthy that, while the bending energy contribution to the elastic energy is small (note the difference in vertical scale in Fig. 3C and Fig. 3D), folding is necessary to achieve a net energy reduction. To demonstrate this, we have plotted the helicoidal elastic energy which corresponds to χ = 1 in Fig. 3C (dashed line). The energy is observed to grow well above the elastic energy of a folded sheet compared to when χ < 1 as in the experiments, showing that folding clearly results in a lower growth in the elastic energy.
We now complement our elasto-geometric analysis with origami construction to explain the folded structure which develops at θ = 180 • . Consider an inextensible sheet (Fig. 4A) which can be folded up or down along the dashed lines, which results in a polygonal spiral origami (Fig. 4B). The apex angle of the isosceles triangular folds is α. An image of an elastic sheet with the same aspect ratio is shown in Fig. 4C, where the thickness of the sheet has been chosen such that it results in the same number of folds as in the origami. We plot the segment angle β from the ȳ-axis made by initial horizontal lines in Fig. 4D. Quantitative agreement is found between the experimental value of β away from the clamps with the expected value (red line) assuming solid body rotation of the triangles, where each fold acts as a hinge. To quantify the role of the strain on the origami pattern, we measure α from shadowgraph images for various strains and find that it decreases with ∆L/L (Fig. 4E). This variation follows from the decrease of λ p and increase in sheet length with strain, if one assumes α = λ p /(L + ∆L). Thus, good agreement can be observed between the origami shape and the twisted sheet away from the clamped edges.
Further quantitative agreement observed with physical cuts mid-way between the clamps can be found in SM 9: Transect Cut Comparisons.
Origami corresponding to spiral accordion folded elastic sheets can be algorithmically generated using α as a parameter. Consider a right angle triangle with height L and angle α/2 (Fig. 5A). (This triangle is also the same as that at the far left side of the sheet in Fig. 4A.) To help understand the geometrical transformation leading to a flat folded origami, we preserve the color when reflecting off a right triangle with respect to its height, and change color when reflecting the triangle with respect to its hypotenuse identified as a fold. A flat-fold origami is thus obtained by applying alternately these two transformations until reaching a given number of folds n. For specific values of α, these transformations result in a regular flat-folded origami where the triangle bases are the edges of a regular polygon with p vertices (Fig. 5A). This polygon is also the convex envelope of a star-shaped polygon composed by the hypotenuses (dashed black line) connecting vertices separated by q consecutive triangle bases (solid orange line).
These origami can be identified by the so-called Schläfli symbols p, q (39), and thus we call them Schläfli origami. By geometric construction, we have p = 180/α and q = (p -1)/2.
Varying the Schläfli symbols (or equivalently reducing the tip angle and increasing the triangle numbers), one can obtain triangle, pentagon, heptagon, and nonagon shaped envelopes. The thickness of the overlapped regions at the center is given by t = (n + 1)t and decreases in integer multiples of t toward the edges.
In Fig. 5B, we show in the first column flat-folded Schläfli origami of higher symmetry, and in the second column their incomplete counterparts obtained by restricting the number of folds to n = 2 W/λ p , keeping p the same. The examples in the 1st, 3rd, and 4th row are symmetric, and the one in the 2nd row is asymmetric, and thus a Schläfli origami with either symmetry can be generated according to our algorithm by starting with a right angle triangle reflecting about the hypotenuse and height equal to the calculated number of folds. Partial Schläfli origami have been denoted with color intensity which is proportional to the number of overlapping domains at the particular location. They can be compared with radiograms of spiral folded elastic sheets which have the same L/W and n (Fig. 5B, third column). The color map in the radiogram is linearly proportional to absorption encountered along the linear path of the x-rays, and thus can be observed to be consistent with those generated by origami.
The correspondence is further quantified by measuring, from the radiograms, the fold angle ψ and α, after a 180 • twist (Fig. 5C-E). α does not vary significantly between the triangles of a given twisted sheet, in accordance with the predictions of the Schläfli origami (Fig. 5E). We find an excellent agreement without any fit parameters for all three measures. Thus, the orientation of the folds is given by our model, which neglects the elastic stretching of the sheet, and reveals the strong connection between twisted sheets, classical Greek geometry and origami.
We now examine the transformation of the folded sheets into yarns by plotting transects at mid-distance between the clamps for θ = 180 We observe that the ratio decreases from 2 to 1 after the secondary instability (denoted by the vertical line), showing a recursive folding of the sheet.
Then, we represent the features by which the multi-layer yarns form by idealized straight crossectional segments connected by curved joints making right (R) or left (L) turns going from one edge of the sheet to the other by the black arrow in Fig. 6A. We take the cross-section shown in Fig. 6A as an example. With this convention, the configuration of the cross-section before and after secondary instability is encoded as LRLR (accordion) and LLRRLLRRL (folded accordion), respectively. The schematics highlights a period doubling by recursive folding in a way which is qualitatively different from twist-less compressed sheets where the sequence of turns LRLR transforms in LRLRLRLR after period doubling (41,42). Fig. 6C show the radiograms at the corresponding θ where the tracked edges of the sheet are marked in red and blue. While the folding and helical wrapping yield complex internal structures, the stretched edges are found to wind around each other, similar to the twisting of two filaments into a rope (43). The crossings between the two edges in the projected plane occur in the yarn-like compact region which starts to develop along the longitudinal axis. To quantify the yarn region, we use the orientation of the segment joining the end-points in the y -z plane to obtain a cross-section orientation angle θ x (Fig. 6D, inset and SM: Fig. S10). We plot θ x in Fig. 6D after the sheet is twisted twice, and find that the twist is localized in a central section L Y , where the local twisting rate θ x /L Y is approximately 4 times greater than θ/L. We use these observations to introduce a geometric model of yarn formation from the accordion folded sheet.
In order to model the growth of the yarns, we assume that the sheet can be divided into three sections with a yarn-like structure of length L Y and two fan-like structures near the clamped edges characterized by fan angle φ, as shown schematically in Fig. 6E. This simplification enables us to retain the fundamental role of the twisted sheet-edge in the elasto-geometric fold model, while circumventing the difficulty in calculating strains.
We measure the evolution of L Y and φ with twist over various L/W , and find that L Y and φ increase with θ (see SM: Fig. S11B andC). The L Y increases quasi-linearly with twist rate depending on the aspect ratio with the onset of yarn formation observed to begin after the secondary instability occurs at θ = θ s . Using trigonometry, we can express the fan angle as a function of the sheet aspect ratio W/L, scaled yarn width D/W , and scaled yarn length L Y /L, yielding:
Before yarn formation L Y = 0, D = 0, and we note from Eq. [5] that tan φ(θ s ) = W/L, where φ(θ s ) can be interpreted as the angle that the diagonal makes with the longitudinal axis. This form predicts an overall decrease of φ with the sheet aspect ratio L/W . We compare φ obtained using various sheets with Eq. [5] in Fig. 6F, and find very good agreement.
Then, the evolution of the yarn length can be understood from the helical wrapping of the fan edges around a cylindrical core of diameter D that encompasses the crosssection of the compacted material in the yarn region. The fan edges are assumed to be in direct contact with the core with an angle φ (see SM 11: Yarn shape analysis and Fig. S11A), thus forming a helix with a local twist rate tan φ/(D/2) where φ remains unchanged upon twist. We further impose the yarn growth rate with twist, dL Y /dθ to be set by the local twist rate, yielding dL Y /dθ = tan φ/(D/2). Using Eq. [5], we find that L Y is modeled by a linear first-order ODE whose solution is:
where Θ = χ (θ -θ s )/(1 -χ). Considering that there are no adjustable parameters, this growth model is in very good agreement with experimental data shown in Fig. 6G.
Thus, we explain the remarkably ordered transformation of flat sheets to scrolled multilayered yarns observed in our experimental study by introducing a series of simplified elasto-geometric models. Since perturbative analysis of the FvK equations reported previously (14,20) cannot be applied to the large deformation regime, we develop a framework combining observed geometry, elasticity, and kinematics to explain the three stages of the transformation.
Beyond the primary transverse wrinkling instability, the Elasto-geometric fold model (introduced in Results Section C) fully explains the torque required to twist the sheet and the observed accordion folding until the sheet comes in self-contact. In this model, the elastic energy is parametrized by decomposing the sheet into longitudinal filaments of length L, and the compaction parameter χ quantifying the folding midway between the clamped edges. Because these two quantities are independent of the choice of the frame of reference, we can compute the elastic energy even when the sheet presents finite rotation, making our model suitable to address finite shape transformation. This is a fundamental difference from FvK models used to explain initial wrinkling instability with twist (20). We find that finite rotation is responsible for the torque saw-tooth profile shown in Fig. 3B, while the precise peak location is modulated by the fold-induced transverse stiffness of the sheet captured by Eq. [4]. In our model, the bending energy stored along the folds can be neglected, an assumption validated experimentally by the excellent agreement between the calculated and measured elastic energy in Fig. 3C.
To explain the flat multilayered structure observed after a 180 • twist, we then introduce a Schläfli origami model in Results Section D, where we consider the inextensible sheet limit (represented as in a sheet of paper) in alternating mountains and valley creases. The resulting folds form a sequence of vertically oriented triangles with complementary orientations. Thus, the observed transformed sheet shape is parametrized by the number of creases and the angle of the triangle apex in our model. In the flat folded state, we show that the origami is twisted by a half-turn forming aregular star-shaped polygons characterized by Schläfli symbols. These flat folded Schläfli origami accurately predict the observed folded structure when the Schläfli symbols and the vertices number are set using sheet L, W , and λ p . Thus, we find the imprint of the primary wrinkling instability persists far beyond onset. Deviation of the shape from our prediction is observed only near the two clamps where significant stretching is needed to satisfy the clamped boundary conditions. Nonetheless, it is remarkable that the kinematics obtained by an inextensible origami model is preserved when significant tension is applied as in the elastic sheets used in our experiments. The key role of sheet stretchability and applied tension in selecting the number of folds, and in organizing the folding into tightly scrolled yarns is thus uncovered by our analysis.
Then, we develop a geometric yarn model in Results Section F to explain the evolution of the folded sheet after a secondary instability. Based on x-ray tomography analysis, we postulate that the structure can be considered as composed of a highly twisted yarn region at the center and weakly twisted fan-like regions connected to the two clamps. Our model is based on simplified kinematics where the edges are straight and co-planar in the fan region, but form two helices winding around each other and a compact cylindrical core idealizing the multi-layered yarn.
Our model explains the decrease in size of two fan-like regions upon twist and the linear growth of the yarn length initially, with an exponential slow-down as the yarn-ends approach the clamps given by Eq. [6]. Thus, we uncover with experiments and complementary framework, the main stages of the transformation of a flat sheet to multilayered yarns.
Finally, it is noteworthy that the PDMS and PVS sheets in our study were used repeatedly because of their hyperelastic nature. This allowed the sheets to be unfolded and reconfigured multiple times during the course of the trials, recording their shapes and torques under different loading conditions. Thus, we demonstrate that our tensional twist-folding strategy can be used to create redeployable functional structures from simple elements with the appropriate choice of materials, an important goal for advanced manufacturing with soft materials (1, 2, 7).
The material properties of the sheets used in the study can be found in the Supplementary Materials (SM) 1: Materials. The shape measurement methods are discussed in SM 2: Optical Imaging, and SM 3: Fluoroscopy and Computed Tomography, respectively.
The data in the analysis corresponding to the measured sheets, twits, and scans can be found in the main document and in Ref. (44) Supplementary Text
The sheet materials used to study twist induced transformations and their properties are listed in Table 1. The PDMS and PVS sheets are hyperelastic, undergoing fully reversible transformations even when stretched by 50% of their physical dimensions. These sheets do not plastically deform and return to their initial flat unstressed state after being twisted and untwisted repeatedly as shown in Fig. S1. The latex sheets show aging effects and were used to perform limited measurements to extend the range of parameters investigated, and examine the robustness of the phenomena. The stress-strain relationship for these materials are similar for the purposes of our study and was shown to be consistent with the Mooney-Rivlin hyperlastic model in previous published reports (44,26).
The sheets were imaged using a digital color PixeLINK PL-D725CU-T camera and illuminated using white light generated with a uniform LED panel. The sheet and the twisting apparatus were placed in between the light source and the camera. Regions with greater thickness absorb and scatter more light, and thus appear darker. In the case of the transparent sheets, refraction of light near the high curvature regions and absorption of light cause the corresponding regions to appear dark. This enables us to locate folds by identifying the change of light intensity across a fold in an image using standard image processing.
A Varian BIR 150/130 Desktop Computed Tomography system is used to noninvasively obtain the x-ray images and scans of the deformed sheets. The three dimensional reconstructed volume arrays of the material density were obtained using 720 radiograms, each consisting of 948 × 688 pixel 16-bit images, in 0.5 degree rotational increments about a central vertical axis. The 3D coordinates of the sheet were then located in each horizontal crosssection to within 100 µm depending on the volume scanned.
We used a common plastic bag made of High Density PolyEthylene with a 10µm thickness. A square of 10 cm wide and 30 cm long is clamped at both ends and stretched by hanging a weight of 4 kg at the bottom rotating clamps resulting in a strain ∆L/L = 2%. The sheet is then twisted by 5 turns (θ = 1800 • ) and held in place while hot air is blown for approximately 2 minutes. The applied temperature is held in the range 80-140 • C which around the melting temperature of the material. When the plastic approaches its melting temperature, it deforms plastically while still offering mechanical resistance. This deformation keeps the fabricated yarn structure intact even after the applied twist and heat are removed.
We introduce a Cartesian coordinate system (x, ȳ, z) attached to the laboratory frame of reference. The material points of the sheet are parametrized by their Cartesian coordinates (x, y) in the load free configuration, projected on the mid-plane (z = 0) (Fig. S2). After transformation, the coordinates of the material points are denoted by (x (x, y), y (x, y), z (x, y)). This parametrization is used to derive elastic quantities (strain, elastic energy, torque, etc). Then, we introduce the parametrization (x, y c (x, s), z c (x, s)) for the shape of the sheet which is better adapted for tomographic cuts. The coordinates of the crosssection at x are denoted by (y c (x, s), z c (x, s)), where s is the curvilinear coordinates. This parametrization is used for the calculation of the mean and Gaussian curvatures. Note that while the components of curvature tensor depend on the choice of the parametrization, the mean and Gaussian curvatures do not.
The geometrical analysis of the sheets are investigated by further processing of the crosssection data obtained with x-ray tomography. We use the morphological operations on binary images available in the Matlab Image Processing Toolbox. The image of a crosssection is first thresholded to obtain a binary image composed of 0s and 1s, and, then, processed by a skeletonized morphological operation using the Matlab bwmorph function. After cleaning up spurious branches, isolated regions and holes, this yields an image with a black background and a 1px thick continuous white line. We construct a list of points M n of the line with a search procedure. The search procedure gives a list that needs to be ordered to reflect the connectivity of a pixel and its neighbor as illustrated in Fig. S4. The information of the connectivity is necessary to compute quantities such as curvatures which are obtained by a spatial differentiation of the neighboring points.
The permutation of the indices to have all points M n ordered according to their connectivity is obtained via the symmetric reverse Cuthill-McKee ordering applied on the adjacent matrix Adj. The adjacent matrix is a sparse matrix where Adj(n, m) = Adj(m, n) = 1 if M n and M m are distinct connected points (i.e. n = m), else Adj(n, m) = 0. If the points are ordered according to their connectivity, the adjacent matrix is tridiagonal. The symmetric reverse Cuthill-McKee ordering provides the permutation triadiagonalizing the adjacent matrix and, thus, provides the permutation to have the points ordered according to their connectivity. Once the points are ordered, we can use finite difference to calculate the quantities such as curvilinear coordinates, gradient, and curvatures.
The surface is parametrized by (x n , y c n (x n , sn ), z c n (x n , sn )), where y c n and z c n are the coordinates of a point M n of a crosssection at xn , and sn is the curvilinear coordinate defined for n > 1 as:
and s1 = 0 (see Fig. S5A). Next, we set an observation window where a small region of the sheet is analyzed. Typically, the portion of the sheet is 10 by 10 pixels, which is a good compromise to have sufficient averaging of the noise while still keeping the spatial variation of curvature.
To calculate the local curvature tensor, it is best to orient the portion of the surface such that the normal is aligned with the z-axis. In this new frame of reference, the shape of the surface can be fitted by a quadratic polynomial of the type:
where κ ij are the components of the curvature tensor. The local normal n is calculated from the moment of inertia I of the surface inside the observation window. Because a surface has a strong geometric anisotropy (the thickness is much smaller than the lateral dimensions), the eigenvector of I with the highest eigenvalue gives the direction of the normal.
The Gaussian and mean curvatures defined as the determinant and half the trace of the curvature tensor are calculated:
and K = det(κ) = κ xx κ ss -κ 2 xs . We can check for the accuracy of the measurement of H using the data corresponding to applied twist θ = 120 • . Fig. S5B shows the crosssection of a PVS sheet with L/W = 3. In this case, the folds are essentially cylindrical, thus the curvature along the fold is much smaller than in the transverse direction. Hence, the mean curvature can be calculated from the radius of the cylinder using ImageJ. The radius of curvature r c ≈ 2.35 mm, thus κ ss ≈ 0.43 mm -1 , and H ≈ 0.21 mm -1 . This is indeed what is obtained looking at the mean curvature map shown in Fig. 2A in the main text.
We develop an elasto-geometric model for the tensional twist-folding of elastic sheets. The model is tractable enough to yield analytical predictions up to self-contact. When self-contacts occur the kinematics becomes complicated but a simplified version of the model, i.e. the Geometric yarn model, can be developed allowing to obtain the evolution of the morphology of the yarns.
In the elasto-geometric model, we neglected the fine morphology of folds and retain only the backbone shape of the deformed sheet. As shown in the main text, the bending energy is much smaller than the stretching energy, at least up to self-contact (see Fig. 3a andb, main document). To obtain the dominant term in the elastic energy, we neglect the bending contribution i.e. B/(Y tW 2 ) → 0.
Because we focus on deriving analytical solutions for the energy and the torque, we restrict ourselves to configurations before self-contact. We consider an elastic sheet of length L, width W , and thickness t stretched along its length by a distance ∆L Young's modulus Y and Poisson ratio ν, as shown in Fig. S6A. Before twist is applied, the initial stretching results in a strain 0 = ∆L/L, and length between clamps L = (1 + 0 )L. To avoid cumbersome expressions, we take ν = 0. With this choice, there is no loss of generality as our findings do not originate from a Poisson ratio effect. As a consequence, the width W and thickness t remain unchanged after loading. The othonormal basis of the Cartesian coordinate system is (x, ŷ, ẑ), where x is oriented along the axis of rotation, ŷ is in the plane of the sheet, and ẑ is oriented along the normal. The coordinate system is centered on the sheet. The material points are parametrized by -L/2 ≤ x ≤ L/2 and -W/2 ≤ y ≤ W/2.
The top clamp initially at x = L/2 is translated by ∆L/2x and rotated by an angle θ/2, and the bottom one initially at x = -L/2 is translated by -∆L/2x and counter-rotated by -θ/2. Then, the top clamp is parametrized by:
and a parametrization of the bottom clamp is obtained by a 180 • rotation with respect to the y axis.
The backbone shape of the deformed sheet (i.e. averaging out the wrinkling and folding structure) is modeled by two ruled surfaces S + and S -, connected at the mid-crosssection (x = 0) along the y axis. As illustrated in Fig. S6A, the surface S + is generated by sweeping transversely a line connecting the top clamp and an effective crosssection at x = 0 and is parametrized by:
where r χ (y) parametrizes the mid crosssection (at x = 0). The parametrization for S -can simply be obtained by a 180 • rotation with respect to the y axis. In the remaining part of this section, we only consider the top half of the surface for our reasoning. Above the onset of the primary instability, the crosssection at x = 0 starts wrinkling and folding. The kinematics is complex and composed of a contraction along ŷ and oscillations in the xz plane, as evidenced in Fig. 3 in the main text. These oscillations are of a much lower extent than the contraction. For the sake of clarity and conciseness, we will neglect these oscillations. We check with numerical simulation that no significant effect on the elastic energy is observed when these oscillations are incorporated in r χ . Then, the effective crosssection at x = 0 is parametrized by:
For a helicoid, the width of the sheet remains constant upon twist, then χ = 1. For a sheet without bending stiffness, the edges of the sheet remain straight. This constraint imposes that χ = cos θ/2. Indeed, by symmetry, the coordinates of the crosssection are obtained as the mean between the top and bottom clamp coordinate, yielding r χ = y cos θ/2ŷ. Fig. S6B shows that the measured χ 1 when the sheet is helicoidal, before the onset of primary instability, and then decreases by an order of magnitude to a plateau as θ → 180 • . We can readily check that the decrease is smaller than cos θ/2 which is explained by the fold-induced resistance to compaction. Measurements indicate that, in absolute value, the rate of change of the compaction parameter is larger near onset of primary instability than near self-contact. In the experiment, the compaction is well described by a decreasing asymmetric sigmoid function:
where Γ(u, a) is the normalized lower incomplete gamma function. This asymmetric sigmoid gives an accurate fit of the data using 4 parameters. θ p is the threshold from primary instability, χ sc the compaction parameter at self-contact. b sets the decrease rate in twist from 1 to χ sc and a sets asymmetry rate of change near θ p and θ sc (twist angle at self-contact). θ p and χ sc are measured. We found a = 4.5, b = 4 adjusting the parameters by the eyes. Using Eqs.[ S1], [S2], and [S3], the parametrization for S + reads:
In the limit where θ is small, we have to first order in θ:
This expression can be compared with the parametrization of a stretched helicoid:
Taking χ = 1, corresponding to a helicoid in Eq. [S6], we have:
Thus, the kinematics introduced in Eq. [S2] coincides with that of a stretched helicoid to first order in θ. Note that Eq. [S8] can be found as the solution of the Föppl-von Kármán (FvK) equations where the deflection is given by r • ẑ = θ L xy. (If the Poisson effect is included, one would have (1 -ν 0 )y ŷ instead of y ŷ. We have simplified the expressions by neglecting the Poisson effect by taking ν = 0.)
Based on the parametrization given in Eq. [S2], the sheet is divided into a series of independent longitudinal strips of width dy which connect the top and bottom clamp with a discontinuous change in slope at x = 0. (In fact, each strip can be thought of as two straight segments connected at x = 0 because we have assumed B → 0.) An individual segment is characterized by the vector u = r c -r χ connected at its two ends and given by:
Further, the length L of the strip is given by L = 2|| u|| and reads:
where, we introduce:
The evolution of F with θ is shown in Fig. S6C using the same compaction parameter as in B and 0 = 0.16. In the limit where F 1, we have:
Further, we have that F ≈ η ≡ θW/L for a slightly twisted sheet (χ ≈ 1 and θ 1). Thus, L = L 1 + η 2 (y/(W/2)) 2 , which is the length of a helix located at a distance y from the axis (20).
The strain along a strip is defined as = (L -L)/L and reads:
In the limit F 1, the strain is given by:
As before, we obtain
for a slightly twisted sheet (χ ≈ 1 and θ 1), which is the strain profile for a slightly stretched helicoid which can be obtained as a solution of the FvK equations. F 2 can be interpreted as a measure of the twist induced strain extended to large morphological transformations.
The tension per unit width acting along the strip is f = Y t û, where û = 2 u/L is a unit vector. The torque M is defined as
The integrand in Eq. S13 is:
Thus we have :
where M 0 = Y W 2 t and ξ = y/W . The torque can be calculated numerically as shown in Fig. S7A. We find that there is a maximum and very good agreement with the experimental data.
We can obtain an analytical expression in the limit where F 1. The torque reads (up to second order in F):
Finally, after integration over the width, the torque is :
From Eq. [S14], we can readily see that the torque is not monotonic due to the sine function. When χ and F do not depend on θ like in the helicoid regime, the location of the torque peak is given by the condition dM/dθ ∼ cos θ/2 = 0, yielding a peak torque at θ = π. For a sheet with no bending-induced stiffness, we have χ = cos θ/2 (see B. Kinematics). M ∼ 2 0 cos θ/2 sin θ/2 = 0 sin θ, where we assume that 0 /12 F 2 /160 which is the regime of high tension explored in our study. In this limit, the torque peaks at θ = π/2. In the intermediate case (1 < χ < cos θ/2), the location of the peak is between π/2 and π which is in agreement with the the torque measurements.
From Eq. [S14], we can also obtain a scaling for the maximum torque admissible by a twisted sheet M max taking χ min ≈ 1/2 and θ p ≈ π as typical values. At high initial pre-strain ( 0 F 2 ), M max ∼ Y W 2 t (W/L) 0 scales linearly with W/L and 0 . In the opposite limit ( 0 F 2 ), M max ∼ Y W 2 t (W/L) 3 is pre-strain independant and grows as the cube of W/L. If the sheet is in the slightly stretched and twisted helicoid regime (θ 1, L ≈ L, and χ ≈ 1), then F ≈ η, the product inside the square brackets also tends to η, and:
Eq. [S15] is the same as in Chopin and Filho (34) using the FvK model. However, within this model, the torque increases monotonically, thus, no maximum is predicted as shown in Fig. S7A.
In our scalar model, the total energy E of the sheet under fixed distance condition is given by:
where E el is the elastic energy and W M the work by the external torque which is given by W M = θ 0 M (θ )dθ , where M is the applied torque.The elastic energy per unit width is given by dE el /dy = 1 2 Y 2 Lt. The stretching energy is then obtained by integration over the sheet width:
It can be noted that this elastic energy corresponds to the membrane energy, as bending rigidity is taken to zero in the zero thickness limit. Further, unlike the FvK equations where configurations are restricted to remain in the small slope limit, finite rotations are taken into account in the model. In the ribbon limit (L/W 1, thus F 1):
with F given in Eq. [S11].
We can provide a complementary definition of the torque based on the elastic energy
where we used the fact that at equilibrium dE dθ = 0. We can check the consistency of the M with Eq.S13. Noting that dE el dθ = dF 2 dθ dE el dF 2 , we have:
and, using Eq. S16,
Inserting Eq. [S18] and Eq. [S19] in Eq. [S17], we recover the expression of the torque in the ribbon regime (L W ) given in Eq. [S14]. Finally, integrating Eq. [S17] with respect to θ, we obtain an expression of the elastic energy in terms of torque:
where
is the stretching energy at zero twist.
The expression of the elastic energy given by Eq. [S20] is independent of the elastic model used. To calculate explicitly E el in our model, we use Eq. [S13]. Whereas, Eq. [S15] can be used to calculate assuming the FvK model. These two calculations can be directly compared with experiments by using the measured torque with a strain gauge. It is worth noting here that the bending modes do not contribute to the torque neither our scalar model nor in the FvK model.
In the FvK thin elastic sheet model, the elastic energy is given by E el = E s + E b , where E s is the stretching energy and E b is the bending energy given by
where B is the bending stiffness, ν the Poisson ratio, H the mean curvature and K the Gaussian curvature. The stretching energy is given by:
The calculated elastic energy E el corresponding to our scalar model, given by Eq. [S16], reduces to E s in the limit where L/W 1, χ ≈ 1, and F ≈ η = θW/L. We start by estimating an upper bound for the bending energy assuming that the sheet folds with a radius of curvature r c . The mean curvature scales as H ∼ r -2 c . The total length of folds is φ f nL, where n is the number of folds and φ f a numerical prefactor correcting for the fact that the folds are slightly tilted and span around 80% of the sheet. The fold surface is therefore 2r c φ f nL where we take 2r c as an estimate for the fold width. Therefore, we have
Substituting B = 1/12Y t 3 , we obtain:
The radius of curvature r c is in the range 2t -D where D is the yarn thickness. The ratio of the estimated bending energy scaled by E 0 (Eq. [S21]) is:
Next, we estimate the stretching energy at a half turn using Eq S14 keeping terms up to second order in W/L. We obtain:
In the case of a helicoid (χ = 1), the elastic energy reduces to:
The change in energy ∆E el = |E el (θ = π) -E hel el (θ = π)| between the helicoid and folded shape is:
We can now express the ratio between the elastic energy estimated in the scalar model using Eq. [S24], and an upper bound of the bending energy estimated in Eq. [S22] using the FvK model as
To have a rough estimate of this ratio, we take typical values for the various parameters (W/t ∼ 10 2 , W/L ∼ 1/2, n ∼ 5, χ ∼ 1/2, φ f ∼ 1, 0 ∼ 0.1, r c /t ∼ 5), and find that ∆E el E b ∼ 10. Thus, folding results mostly in a reduction of the stretching energy associated with a modest increase in the bending energy, typically an order of magnitude smaller. This rough estimate is consistent with the results shown in Fig. 2H and I in the main text. At θ = 270 • , ∆E el ∼ 0.5 E 0 while the bending energy is E b ∼ 0.05 E 0
To further test the correspondence between the origami model and the twisted sheet by a halfturn, we consider a physical cut midway across the folded sheet, which results in a sawtooth shape which is related to the number of folds (Fig. S8F). In fact, it is possible to derive a recursive relation for the cut angle γ m (Fig. S8G) which can be related to the triangle number m from the center as
where, α = λ p /L. Then, by measuring the average slopes between each peak and valley, we compare the ones observed in the folded hyperelastic sheet, the corresponding origami, and prediction of Eq. S26 in Fig. S8C. Good correspondence can be observed with systematic lower slopes because of relieving of the stretching due to cutting of the sheet.
We characterize the highly deformed shape of a twisted sheet when a scrolled yarn has formed. As illustrated by the fluoroscope shown in Fig. S9A, a projected sheet crosssection can be observed to vary considerably along its length. Hence, we measure the degree of compaction of the sheet as a function of x by calculating the diameter D x of an enclosing circle by considering each scanned crosssection as illustrated in Fig. S9B. The center of the circle which envelops the sheet is located at the center of mass of the sheet crosssection, and the circle radius D x /2 is given by the distance from the center to the farthest point in that crosssection. From the variation of D x observed in Fig. S9C, we conclude that the sheet can be considered to be composed of a central compact yarn-like region and two triangular fan-like regions with bases located at the respective clamps. We first discuss two bounding values of D x based on the observed variation. A helicoidally deformed sheet corresponds to the upper bound, D = W , and χ = D/W = 1. To estimate the lower bound, we note that the initial area of the crosssection is tW . Assuming a circular yarn without holes, the area is π 4 D 2 . Assuming crosssectional volume conservation in each slice, tW = π 4 D 2 . Thus, the minimal compaction parameter is χ = 4t πW . For the sheet considered here, then minimal χ = 0.083. In the central region, we measure χ = 0.11, which is reasonable considering that we have neglected out-of-plane motion of the sheet in each crosssection as well as Poisson effect, and the fact that the sheet may not be fully packed.
After examining the envelop of the sheet, we now investigate the internal structure by tracking the two unclamped edges as a function of applied twist for θ = 180 • , 360 • , 540 • , and 720 • . We locate the coordinates of the two edges (indicated by a blue and red markers) in each yzplane with image analysis software ImageJ (Fig. S10A). Then, we calculate the distance D x and orientation angle θ x of the segment joining those two points as a function of x. In Fig. S10B and C, we plot the distance D x and the angle θ x . Although the minimal distance between the two edges min(D x /W ) = 0.04 is smaller than the corresponding χ, we find the same trend looking at D x as compaction parameter D x shown in Fig. S10C. This is to be expected as the two edges can actually be in contact. Thus, according to parametrization, the sheet develops two distinct regions with increasing twist, in which the central compact yarn-like region grows at the expense of the fan-like regions near the two clamps.
We model the scrolled yarn formation as a structure with a core of width χW around which the fan edges wrap around (see Fig. S11A). The fan edges are assumed to be in direct contact with the core with an angle φ which remains unchanged upon twist. Thus, the fan edges transition into a helical structure wrapped around a cylinder with the same angle φ in our model. Since φ evolves with θ, the helix has a pitch which continuously evolves with θ.
To obtain the growth mechanism, we write that L Y increases at a rate given by τ Y which is the twisting rate measured at the yarn ends. Inside the central region, the twisting rate is on average ∆θ/L Y but varies along the longitudinal direction as may be noted from Fig. S9C. The rate of yarn growth is then given by: dL
For a constant twist rate τ Y , we obtain a linear increase of the yarn length with twist, L Y = ∆θ/τ Y , where ∆θ = θ -C, and C is a constant of integration. However, as the yarn increased in length, we observed that the fan angle also increased. Thus, we expect the twisting rate -which measures the rate at which the sheet wraps around the yarn per unit length -to increase as well. This increase is suggested by the slope variation in Fig. S9C. Focusing on the region near the yarn ends, we can express τ Y in terms of the fan angle φ and compaction parameter χ as (Fig. S11A):
This expression captures the variation of the twist rate as the yarn grows. We can see that as the fan angle increases (as observed experimentally), the twist rate is increasing as well.
Then we obtain an expression of the yarn length growth as a function of φ and D, dL Y /dθ = (D/2)/ tan φ. Using the expression of φ found previously, we have:
Thus, we obtain a first order linear ODE for the evolution of the yarn length
.
where L Y (θ s ) = 0, i.e. the yarn starts growing after the second instability occurs. The solution of the equation is:
As shown in Fig. 6 in the main document, when the yarn length is rescaled by L and plotted against Θ = χ (1-χ) (θ -θ s ), we observed a very good collapse of the data on the predicted master curve given by Eq. [S27]. The material coordinates of the mid-section (z = 0) in the load free configuration are (x, y) and coincide with the laboratory coordinates (A). After transformation (B), the material coordinates are denoted (x (x, y), y (x, y), z (x, y)), and are functions of the material coordinates in the load free configuration. Tables S2.
Materials