Frustration Propagation in Tubular Foldable Mechanisms

Shell mechanisms are patterned surface-like structures with compliant deformation modes that allow them to change shape drastically. Examples include many origami and kirigami tessellations as well as other periodic truss mechanisms. The deployment paths of a shell mechanism are greatly constrained by the inextensibility of the constitutive material locally, and by the compatibility requirements of surface geometry globally. With notable exceptions (e.g., Miura-ori), the deployment of a shell mechanism often couples in-plane stretching and out-of-plane bending. Here, we investigate the repercussions of this kinematic coupling in the presence of geometric confinement, specifically in tubular states. We demonstrate that the confinement in the hoop direction leads to a frustration that propagates axially as if by buckling. We fully characterize this phenomenon in terms of amplitude, wavelength, and mode shape, in the asymptotic regime where the size of the unit cell of the mechanism~$r$ is small compared to the typical radius of curvature~$\rho$. In particular, we conclude that the amplitude and wavelength of the frustration are of order $\sqrt{r/\rho}$ and that the mode shape is an elastica solution. Derivations are carried out for a particular pyramidal truss mechanism. Findings are supported by numerical solutions of the exact kinematics.


I. INTRODUCTION
Origami tessellations are sheets of paper folded following a repetitive crease pattern.Beyond their artistic value, they now constitute a prototyping platform for a remarkable array of foldable and deployable structures useful in engineering [1].Origami tessellations are often modeled as linkages: assemblies of rigid elements representing the facets of the origami, hinged along edges representing the crease lines.From this point of view, the kinematics of origami are no different than for truss mechanisms.Thus, more generally, we will call tessellation any linkage that is periodic (i.e., that is invariant by translation along a two-dimensional lattice of vectors) and we will call folding any inextensional deformation that is compatible with the linkage kinematics.Suppose now that the tessellation undergoes a uniform folding, i.e., a folding where fold angles remain periodic even if the linkage does not.Then, generically, the tessellation embraces a curved, often cylindrical, midsurface.See the examples of the waterbomb tessellation, the Yoshimura pattern and the Ron-Resch pattern [2] as well as the pyramidal truss mechanism illustrated on Figure 1a.Hence, it appears that the in-plane stretch, or contraction, due to the folding of a single unit cell, is coupled to an out-of-plane bending motion that is necessary to maintain geometric compatibility among consecutive unit cells.We refer to tessellations that conform to this description as coupled.Exceptions exist however: there are tessellations that fold uniformly while maintaining a planar midsurface; we refer to them as planar.Examples of planar tessellations include the Miura-ori (Figure 1b), the eggbox pattern and Chebyshev nets.
Planar tessellations have received considerable attention in recent years and have been studied in the context of "geometric mechanics" by several authors.Schenk and Guest [6] and Wei et al. [7] computed the Poisson's coefficient of the Miura-ori and showed that it is equal to the ratio of normal curvatures observed in its anticlastic bending.Similar results were then obtained for the eggbox pattern [3], for the "morph" pattern [8], and, more recently, for a whole class of smooth and polyhedral surfaces of translation with straight or curved creases [5,[9][10][11].Nassar et al. systematically leveraged the Poisson's coefficient identity, and classical differential geometric tools, to compute the tubular states folded out of the Miura-ori [4], the eggbox pattern [3,12] and the "Mars" tessellation [5].Other origami tubes were designed and investigated by Tachi and Miura [13] and by Filipov et al. [14] with focus on axial deployment.However remarkable, the foregoing results can hardly be extended to coupled tessellations mainly because of one difficulty: size dependence.Indeed, the midsurface geometry of a coupled tessellation greatly depends on the size of the creases: finer crease patterns produce tighter midsurfaces, with divergent curvatures; see, e.g., Figure 1c.Similarly, the waterbomb tessellation, the Ron Resch pattern and the Yoshimura pattern all admit pairs of flat-unfolded and flat-folded states, but the paths from one to the other necessarily go through curved states whose curvatures are inversely proportional to the size of the creases.By contrast, the midsurface geometry of a planar tessellation, say the Miura-ori, depends far more on the folding angles than on the size of the creases: finer crease patterns rapidly converge to a limit, non-singularly curved, surface; see, e.g., Figure 1d.[3][4][5]).Note that (c) and (d) feature non-uniformly folded states, i.e., states where folding angles vary in space.
The purpose of the present paper is to alleviate to some extent this difficulty: we propose a differential geometric framework for the characterization of the curved midsurfaces embraced by coupled tessellations under non-uniform foldings, and provide appropriate asymptotic scalings for their size-dependent geometry.The framework is relevant for tessellations where the size r of the unit cell is small compared to the typical radius of curvature ρ of the embraced surface.Analysis of the folding of tubular states in particular reveals that confinement in the hoop direction causes bulges to develop at equal intervals in the axial direction, a phenomenon here referred to as "frustration propagation".This phenomenon has recently been studied numerically by Imada and Tachi [15] as a dynamical system with an area-preserving quality that explains periodicity.Here, focus is on the size-dependence of the frustration.Specifically, we show that, generically, the amplitude and wavelength of the frustration are of order r/ρ and that the profile of the frustration is an elastica solution.Earlier hints of the size dependent behavior of coupled tessellations can also be inferred from the boundary layers observed in [3] and in [2].
The paper begins by introducing one particular coupled tessellation in the form of a pyramidal truss mechanism.Uniform foldings and corresponding stretch-bend kinematics are explored first.Then, non-uniform foldings are described using the differential geometry of the midsurface.The key element is an equation that determines the metric of the midsurface in function of its curvatures and of a characteristic length scale.The theoretical implications are explored for tubular states and the main results regarding frustration propagation are proven.Numerical solutions of the exact kinematics are shown to match the theoretical predictions.Consider the pyramidal truss mechanism of Figure 1a initially introduced in [3].Its unit cell is a spherical four-bar linkage with a single Degree Of Freedom (DOF) assigned to either of the two central angles, θ and θ * ; see Figure 2. Let O designate the apex of a pyramid with side length r and let (A, B, C, D) be the vertices at its base.Let (e 1 , e 2 , n) be an orthonormal basis such that e 1 and e 2 are respectively aligned with # » DB and # » AC.Then, with O as the origin, the remaining vertices are given by the coordinates with Last, θ and θ * depend on one another through the constraint AB = 1 which yields A uniform folding is characterized by three elements: angle θ and two rotations L and R that map a unit cell to two of its neighbors.Compatibility requires that L and R commute as shown on Figure 3a.Therefore, L and R share the same axis of rotation given by some unit vector t.Furthermore, Thus, compatible axes of rotation t must be in the bisector plane of the diad ( # » AB, # » DC) as well as in that of ( # » AD, # » BC).These two planes contain e 1 and e 2 and are, in fact, the same.Therefore, the axis of rotation is given by some angle φ such that t = e 1 cos φ + e 2 sin φ. ( Once φ is given, R and L are uniquely determined through (4).By iterating L and R, the whole folded state can be constructed out of a single unit cell.Given that L and R share the same axis t, and that t is invariant under the action of L and R, it comes that the folded state is a cylinder of axis t.In conclusion, uniformly folded states constitute a 2-DOF family of cylinders parametrized by θ and φ: the former prescribes the folding angles within a unit cell whereas the latter prescribes the folding angles between neighboring cells.For other tessellations that embrace similar cylindrical states, see [2,16].
To gain further insight, let us explore the particular case φ = 0, i.e., where the axis of rotation t aligns with e 1 .Then, by symmetry, whatever rotation maps In other words, R = L.The common angle of rotation denoted α can be deduced algebraically from equation ( 4) or geometrically from Figure 3(b, c), namely Hence, from the same Figure, it comes that the radius of the embraced cylinder, as measured from the cylinder axis to the base of a pyramid, is More generally, consider how a unit cell paves a cylinder by the repeated action of rotations as shown on Figure 4a.By inspection of Figure 4b, it comes that the radius of curvature is where α is the angle of rotation R and is given by see Figure 4c.Expanding the various dot products leads to A generic uniform folding is illustrated on Figure 5a.Profiles of the normalized curvature r/ρ v.s.angle φ are depicted on Figure 5b for various folding angles θ.The curvature is smallest (globally or locally) at φ = 0 and at φ = π/2; that is: the least curved states are cylinders of axes e 1 or e 2 .Furthermore, the radius of curvature ρ becomes comparable to r as soon as θ departs significantly from π/2.For θ = π/2, there is a singularity: in that case, all φ ̸ = ±π/4 lead to the same planar state with zero curvature.As for φ = ±π/4, the axis of the cylinder aligns with either side of the base and radius ρ can take any value larger than r.

B. Linearization
The derivations of the above section show that any extension in direction e 1 or e 2 , as measured by angle θ, is coupled to bending about some axis t with a radius of curvature ρ.Furthermore, ρ scales like r.In other words, for equal θ, finer tessellations wind more tightly.The divergence of the curvature in 1/r towards infinity as r decreases can only be avoided by limiting the magnitude of the folding angle θ to values close to π/2.Thus, henceforth, focus is on states such that where γ represents a small relative extension.The idea behind this restriction is that it enforces r ≪ ρ in a way that makes it possible to speak of a smooth midsurface for more general, non-uniformly folded, states.However, it is not clear how small γ should or will be compared to r/ρ and both are kept as small independent parameters for the time being.Then, to leading order, the radius of curvature is Beyond ρ, it will prove convenient to compute how the midsurface bends relative to directions e 1 and e 2 .Specifically, let Therein, E is the normal change in e 1 transported to the next cell over in direction e 1 across a distance DB; thus, it quantifies the normal curvature of the midsurface in direction e 1 ; refer to Figure 3a for notations.Similarly, G quantifies the normal curvature in direction e 2 and F quantifies the torsion of the midsurface.The angles of rotations L and R, called α and β respectively, are given by equation ( 4) and read to leading order in γ.Therefore, It appears then that (Er/γ, Gr/γ, F r/γ) depend solely on φ meaning that these quantities must satisfy two algebraic constraints, namely These algebraic constraints define the accessible bent states as prescribed solely by the angular extension γ, regardless of φ ̸ = ±π/4.The latter constraint, in particular, says that the Gaussian curvature vanishes, which is expected since cylinders have zero Gaussian curvature.The case φ = ±π/4 has been avoided above but, should it arise, the restriction r ≪ ρ implies γ = 0 and both algebraic constraints remain valid.Now since these states all reside in the vicinity of the planar state (i.e., θ = π/2) and are accessible, a linear superposition should also be accessible, even if non-uniformly folded.For instance, by combining a state A bent about e 1 (i.e., φ A = 0 and γ A ̸ = 0) with a state B bent about e 2 (i.e., φ B = π/2 and γ B ̸ = 0), a doubly curved state is obtained that is extended through γ = γ A + γ B and whose curvatures are given by E = E B , G = G A and F = 0. Evidently, such a combined state will no longer be cylindrical or even have zero Gaussian curvature (i.e., EG − F 2 ̸ = 0) but, remarkably, it will still satisfy the first constraint, namely G − E = 2γ √ 2/r.More generally, this constraint holds for uniformly folded states and is linear, therefore it holds for any linear combination of folded states.By contrast, the second constraint, namely EG − F 2 = 0 holds for uniformly folded states but does not carry over to more general states.

C. Differential geometry of the midsurface
We are now ready to make the scale transition in kinematic modeling from the discrete level to the continuum level.Let (ξ 1 , ξ 2 ) → x(ξ 1 , ξ 2 ) parametrize the midsurface of a general, non-uniformly folded, state.The in-plane deformations of the midsurface can be quantified using the metric tensor g = g ij with Here, we identify which implies, to first order in γ, that Above, and in the following, it is understood that γ ≡ γ(ξ 1 , ξ 2 ) can depend on space coordinates so as to allow for the continuum description to encompass non-uniformly folded states.The only restriction in that regard is that variations in γ should occur over length scales larger than r.
As for out-of-plane deformations, they are captured by the second fundamental form b = b ij where Here too, we identify b 11 ≡ E, b 22 ≡ G and b 12 ≡ F so that the constraint is systematically enforced point-wise, i.e., at all positions (ξ 1 , ξ 2 ).Again, the other constraint, namely det b = 0, only holds for uniformly folded states (i.e., γ = cste) but not in general.Solving for the midsurface then amounts to finding surfaces whose metric and second form, g and b, satisfy equations ( 19) and ( 21).Alternatively, it is possible to combine both equations into one statement: the tessellation embraces midsurfaces whose metric depend on the second form through It is worthwhile to highlight that for several planar tessellations, including the Miura-ori, the eggbox pattern as well as other patterns with no stretch-to-bend coupling, the coefficients of the second fundamental form b still satisfy a linear metric-dependent constraint similar to (21) albeit one that is homogeneous, e.g., of the form p(γ)b 22 − q(γ)b 11 = 0; see, e.g., [3][4][5].Here, by contrast, the foregoing constraint exhibits a "source" term in γ/r that is both metric-and size-dependent.

A. Geometry of tubular states
We define a tubular state as one that embraces a midsurface of revolution.Specifically, let where ρ and z are to-be-determined functions, q ≡ 2π/W and W is the width of the tessellation in the flat reference state; see Figure 6a.Then, computing x 1 , x 2 and so on, it comes that where a prime denotes d/dξ 2 .The key equation ( 22) then yields two Ordinary Differential Equations (ODEs), namely In the following, we solve these equations in the asymptotic regime that is consistent with the adopted continuum description, i.e., in the limit where γ → 0 and r/ρ → 0.

B. Asymptotics of size dependence
We deal with equations (25) as algebraic equations to begin with.Then, up to an error smaller than γ, and there exists a function ω ≡ ω(ξ 2 ) such that since, now to leading order in γ, we have ρ ′2 + z ′2 = 1.Substituting back into expression (24) for γ then yields where terms in γ 2 and in γqr = O(γr/ρ) have been neglected, both being small relative to γ. Taking the derivative d/dξ 2 and replacing γ ′ with its leading-order expression in terms of qρ ′ leads to But ω ′ is of order b 22 ≪ 1/r meaning that, to leading order, ω is solution to the second-order ODE Therefore, it appears that ω oscillates at a length scale Accordingly, it is appropriate to re-scale the ξ 2 coordinate and introduce a function ω such that Then, ω is solution to the elastica ODE, also known as the simple pendulum ODE, ω′′ + sin ω = 0.
Note that further linearization of the elastica ODE is not warranted in general since there is no reason for ω to be small unless further assumptions regarding the smallness of the initial conditions ωo and ω′ o are made.In any case, once a solution to the elastica ODE is chosen by setting ωo and ω′ o at, say, ξ 2 = 0, the profile of the state can be determined by integration as in or, by re-scaling, By analogy with the simple pendulum (Figure 6b), we know that sin ω, a quantity analogous to the horizontal deflection of the pendulum, oscillates periodically with zero average.Hence, ρ oscillates periodically with an amplitude of order L = O( √ rW ).By contrast, cos ω, a quantity analogous to the vertical deflection of the pendulum, oscillates periodically but does not necessarily have zero average.Hence, generically, z grows linearly.In conclusion, tubular states exhibit a periodic frustration whose wavelength and amplitude are of order √ rW and whose profile is an elastica solution.
Note that the initial value z o is arbitrary and inconsequential because it amounts to a rigid body translation along the axis of revolution.As for ρ o , it can be obtained by reconsidering the expression of γ, namely where terms of order qr = O(r/W ) have been neglected in favor of terms of order O( r/W ).This provides an alternative expression for ρ, that is and with it an expression for ρ o in function of ω′ o .Note that, for the sake of consistency, the initial condition ω′ o must be chosen of order O(1), i.e., small relative to L/r = O( W/r).
Last, it is worthwhile to shed some light on the normal curvatures b 11 and b 22 .On one hand, ω = arctan(ρ ′ /z ′ ) is the incidence angle of the tangent plane relative to the axis of revolution.Its derivative can directly be interpreted as the normal curvature in the axial direction, specifically This normal curvature blows up in the limit r → 0 and is of order 1/L = O(1/ √ rW ).On the other hand, the normal curvature in the hoop direction is and remains finite, being of order 1/W .

C. Phase diagram
The elastica ODE admits an invariant analogous to the mechanical energy of the simple pendulum The iso-contours of I define the solution orbits in (ω, ω′ )-space, each orbit being in correspondence with a tubular state; see Figure 6c.Two orbits in particular, namely I = ±1, degenerate into two fixed points corresponding to (ω = 0, γ = 0) and (ω = π, γ = 0); these are two cylindrical states with the pyramids pointing outward and inward, respectively (Figure 6d, h).Since γ = 0 in both cases, it is insightful to seek the next order in its asymptotic expansion.Going back to equation ( 24) and letting ρ be a constant lead to where the sign is positive for ω = 0 and negative for ω = π.In other words, achieving a cylindrical state requires closing the pyramids slightly (as seen in the axial direction) when the pyramids are pointing outward and opening them slightly when they are pointing inward.Interestingly, these two fixed-points have different stability properties.The cylindrical state with the pyramids pointing out is stable: deviations lead to oscillations that grow throughout the profile at once (Figure 6d-f).By contrast, the cylindrical state with the pyramids pointing in is unstable: deviations lead to oscillations growing at localized intervals that, should self-penetration be precluded, cannot exist except at the edges of the tubular state (Figure 6g-i).These edge states correspond to a pendulum that remains near the top (unstable) equilibrium position most of the time and that, however swiftly, swings by the bottom (stable) equilibrium position.

D. Numerical solutions and error analysis
The tessellations produced throughout the paper and in Figure 6 in particular are rendered by solving the exact discrete kinematics.The solution algorithm proceeds iteratively by constructing the nodes of the folded state one "zigzag" at a time as illustrated on Figure 7a.Mainly, given the nodes that belong to zigzag number i, the nodes of zigzag number i + 1 can be constructed by completing the bases of the pyramids that are enclosed between zigzags i and i + 1, one pyramid at a time.This is possible because each enclosed pyramid has a unique degree of freedom.Further detail can be found in [3,4].The algorithm is implemented in Python and its code is available online (see Supplemental Data).Here, rotation symmetry is leveraged and the algorithm can be initialized simply by giving the three parameters (r/W, θ o , ω o ), where θ o = π/2 − γ o .Profiles of γ are computed for small but finite r by solving the discrete kinematics and are shown on Figure 7b.The profile of γ obtained by solving the elastica ODE is shown as well.The plots demonstrate that the discrete solution, hereafter denoted γ discrete , converges to the continuous one, hereafter denoted γ ODE .The convergence error is further computed over a normalized r-independent range 0 ≤ ξ 2 /L ≤ T and is plotted on Figure 7c.Note that the γ profiles have been normalized as well as anticipated by equation (36).The logarithmic scale shows that the error e decays with the side length r like r 1/2 .

IV. CONCLUSION
The paper proposes a differential geometric framework that allows to deal with the size-dependent kinematics of coupled tessellations, i.e., of periodic foldable structures where in-plane folding is coupled to out-of-plane bending.The framework succeeds in predicting size effects both qualitatively and quantitatively for tubular states where hoop confinement produces a frustration that propagates axially.The convergence of the framework is somewhat slow with the error being of order r 1/2 where r is a small size parameter.Thus, improvements that take higher order corrections into account are desirable as they can improve convergence speed.In principle, such corrections will intervene at two places: in the linearization of the discrete kinematics and in the asymptotics of the governing non-linear ODE.
The main conclusion of the analysis, namely that axisymmetric frustrations are elastica-shaped periodic with wavelength and amplitude of order r 1/2 , is not specific to the current pyramidal truss mechanism and should be typical of coupled tessellations that satisfy two properties: (i) the coupling is linear, meaning that g = g(b) is linear (affine strictly speaking); (ii) the tessellation is rectangular so that the coordinate lines aligned with the axes of symmetry deform without shearing.Under these conditions, equation ( 19) is generic and is sufficient for our conclusion to hold.That being said, extensions to other potentially non-linearly-coupled or oblique tessellations are of interest as well and are yet to be investigated.

FIG. 2 :
FIG.2: Annotated unit cell of the pyramidal truss mechanism: the truss is equilateral with side length r.

)FIG. 3 :FIG. 4 :
FIG. 3: Analysis of uniform foldings.(a) Compatibility of rotations: arrows denote transitions between neighboring unit cells and not the axes of rotation.(b) Axis and angle of rotation in the particular case φ = 0. (c) Another view of (b) in the plane normal to t showing the radius of the cylindrical state: reported angles and lengths are measured in the normal plane assuming r = 1.For the full cylindrical state, refer to Figure 1b.

FIG. 6 :
FIG. 6: Tubular states.(a) Notations.(b) Simple pendulum analogy.(c) Phase diagram: iso-contours of the analogue to mechanical energy I. (d-i) Two views of various tubular states for increasing I: (d, h) are uniformly folded states; the other states are non-uniformly folded.

FIG. 7 :
FIG. 7: Convergence analysis.(a) Schematics of the discrete solution algorithm: knowledge of zigzag i suffices to uniquely determine zigzag i + 1.(b) Normalized folding angle γ v.s.normalized discrete curvilinear coordinate as measured by a nodal index i: curves are obtained for decreasing r (see legend); the limit r → 0 is the solution to the elastica ODE.(c) Convergence speed: Data points show the mean quadratic error e as obtained from (b); the seamless line shows the trend e = r 1/2 .Analysis is carried for ωo = π/2, ω′ o = 0. Other states show similar trends.