Simulation and design of isostatic thick origami structures

Thick origami structures are considered here as assemblies of polygonal panels hinged to each other along their edges according to a corresponding origami crease pattern. The determination of the internal actions caused by external loads in such structures is not an easy task, owing to their high degree of static indeterminacy, and the likelihood of unwanted self-balanced internal actions induced by manufacturing imperfections. Here we present a method for reducing the degree of static indeterminacy which can be applied to several thick origami structures to make them isostatic. The method utilizes sliding hinges, which permit also the relative translation along the hinge axis, to replace conventional hinges. After giving the analytical description of both types of hinges and describing a rigid folding simulation procedure based on the integration of the exponential map, we present the static analysis of a series of noteworthy examples based on the Miura-ori pattern, the Yoshimura pattern, and the Kresling pattern. The method can be applied for the design and realization of thick origami structures with adequate strength to resist external actions.

In contrast to thin origami structures, thick origami structures are assembled from a certain number of polygonal panels, hinge-connected to each other along their edges, according to a given crease pattern.Thick origami structures are typically overconstrained, that is, they are statically underdetermined structures, or, in other words, they possess several self-stress states.This poses the problem of the likelihood of unwanted self-balanced internal actions induced by manufacturing imperfections.At the same time, the choice of constitutive relations for the internal actions associated with well defined strain measures is not trivial.For these reasons, to perform the structural design of these systems is a challenging task.
In this work, we detail a strategy for reducing the degree of static indeterminacy of thick origami structures, possibly bringing it down to zero, in order to make them isostatic structures.We take advantage of the introduction of sliding hinges, in addition to conventional door hinges, in a rigid origami model.Sliding hinges provide an additional degree of freedom with respect to door hinges, in that they permit the relative translation between adjacent panels along the shared hinge axis, other than the relative rotation between the two panels about the same axis.Therefore, by replacing a door hinge by a sliding hinge in a thick origami realization, either one degree of static indeterminacy is removed or one degree of kinematic indeterminacy is added.We found several noteworty cases in which the present strategy is successful in making the thick origami structure isostatic, so that the internal actions induced by external loads can be uniquely computed.
Preliminary result of our procedure were presented in [25] for the case of a Yoshimura crease pattern.Here, other than providing the details of the modeling equations and the simulation algorithm, we present results concerning the statics of thick Miura-ori, Yoshimura, and Kresling origami.
In the following, after recalling useful counting rules for the degrees of kinemaic and static indeterminacy (Section 2), we give an exact finite kinematic description of the motion of rigid origami structures with door hinges and sliding hinges, together with the infinitesimal counterpart (Sections 3.1 -3.3).Folding/deployment simulations can then be performed by numerical integration of the exponential map, by the procedure outlined in Section 3.4.The equilibrium equations relating external loads and internal actions are obtained by duality (Section 4.1).Afterward, the equilibrium equations are solved numerically for selected examples to demonstrate the usefulness of the method (Section 4.2).We close by discussing the obtained results and future extensions (Section 5).

Counting mechanisms and self-stress states
In this section, some handy results on the number of internal mechanisms and selfstress states of rigid origami structures, which corresponds respectively to the degrees of kinematic and static indeterminacy, are reviewed for later use (cf.[25]).
A starting point is provided by the convex triangulated polyhedron theorem [26].The theorem states that a pin-jointed bar framework constructed on a convex triangulated polyhedron, by placing the bars and pins of the framework on the edges and vertices of the polyhedra, is isostatic.Then, by removing one edge from such a framework, an internal mechanism is produced, while the four edges surrounding the removed edge form the boundary of a so-called hole [27].Next, by removing one of those four edges, a bigger five-edge hole is formed and another internal mechanism is introduced.Hence, the number of mechanisms introduced by a hole bounded by E b edges is E b − 3.By observing that a origami with a triangulated crease pattern has the same edge-vertex connectivity of a convex triangulated polyhedron possessing one hole with a number of boundary edges equals to the number of boundary edges of the origami, one can conclude that a pin-jointed framework constructed on a triangulated origami crease pattern, in a non-singular configuration, has E b − 3 independent internal mechanism.
In order to determine the number of self-stress states, one can consider again a framework built on a convex triangulated polyhedron and add one edge between the two distant vertices of two adjacent triangles.In this case, the four vertices of the two adjacent triangles form a so-called block, and the framework acquires one self-stress state.By matching the polygonal panels of the origami with the blocks of the pinconnected framework, it is easy to see that each polygonal panel of the origami with V vertices corresponds to V − 3 independent self-stress states in the framework.
By considering the kinematics of a panel-hinge model, in which only the relative rotation about the common edge between adjacent panels is permitted, one can arrive at the same conclusions of the bar framework model.However, the predictions regarding the statics differ in the two models.In fact, it is easy to see that the panel-hinge model gives a higher number of self-stress states.In particular, it was shown in [25] that S P H = S BF + 3 V i , with S P H and S BF the degree of static indeterminacy computed according to the panel-hinge and bar-framework models, respectively, and V i the number of the internal vertices of the origami.Figure 1 illustrates these quantities for the square twist origami.

Kinematics
In this section, we first derive the exact constraint equations for sliding hinges and door hinges between adjacent panel-shaped bodies, together with the corresponding linearized versions.Then we give them an equivalent description by composition of point slider and spherical hinge constraints.After that we present an integration procedure of the kinematic equations based on the Newton-Raphson algorithm.In the following, we consider the motion of a body B in the Euclidean space with respect to a fixed reference frame as described by the rigid motion g = (R, u) in terms of a rotation tensor R ∈ Orth + and a translation vector u ∈ V .
Let B 1 and B 2 be two bodies connected to each other by a hinge with axis r in the reference configuration, and let r 1 and r 2 the images of r in the rigid motions g 1 = (R 1 , u 1 ) and g 2 = (R 2 , u 2 ) of B 1 and B 2 respectively: with p is the position vector of a point P on r, and t is the unit vector parallel to r.

Sliding hinges
In case of a sliding hinge, we have that r 1 = r 2 , that is, that is with s 1 → s 2 (s 1 ) and s 2 → s 1 (s 2 ) and s 1 • s 2 the identity.By differentiating with respect to the parameter s 1 the first of these two equivalent relations, we get s ′ 2 R 2 t − R 1 t = 0 , so that s 2 ′ = ±1 and R 2 t = ±R 1 t.Assuming that at time t 0 , R 2 (t 0 )t = R 1 (t 0 )t, the continuity of the motion implies that with α an arbitrary vector.Then, substitution in (1) yields α × (p 2 − p 1 ) = 0 , where p 1 and p 2 denote the image of the position vector p under the rigid motios g 1 and g 2 of the two bodies: The constraints (1) can be realized, for example, by the scalar equations or, by the scalar equations, where n, m are two linearly independent unit vectors orthogonal to t.
We choose as constraint equations the arithmetic mean of equations of the two groups.Consequently, the constraint functions are The time derivatives of the constraints functions have the expressions To obtain these expressions we have used the differentiation rules with W i = Ṙi R T i and w i = ui − W i u i .By introducing the angular velocities ω i as the axial vectors of W i , i = 1, 2, (2) are expressed as follows Then, the first order approximation of the constraint system in the neighborhood of the configuration (R 1 , u 1 , R 2 , u 2 ), has the expression

Door hinges
For door hinges we start with the condition With the notation used in the previous section, t 2 − t 1 = 0 (as in the case of a sliding hinge) , Then the constraint function are The time derivatives of the constraints functions have the expressions or, in terms of angular velocities, Then, the first order approximation of the the constraint system in the neighborhood of the configuration (R 1 , u 1 , R 2 , u 2 ), has the expression with * p i the axial tensors of p i , i = 1, 2, and I the identity tensor.

Point constraints
We express here the hinge constraint between two bodies in terms of point constraints, that is, a sliding hinge is obtained with two sliders located on the hinge axis, a door hinge is obtained as a fixed point and a slider located on the hinge axis.

Single point slider
Let the axis of the guide be described in the reference configuration by the line r solidal to B 1 that passes through the points O + p and O + q ; and let O + q be the point of B 2 constrained to slide on the guide.Under the motion of the two bodies the vector positions p, q are transported to the vectors and the line r to the lines with t = vers(q − p) and i = 1, 2.
The constraint imposes that the point p 2 be the position vector of a point of r 1 , that is It follows that the constraint equations are as follows: The derivatives with respect to time of these functions have the expressions Then, the first order approximation of the constraint equations in the neighborhood of the configuration (R 1 , u 1 , R 2 , u 2 ), has the expression , and

Sliding hinge as a double point slider
With the notation introduced previously, the constraint imposes that the points p 2 and q 1 be the position vector of a point of r 1 and r 2 , respectively; then It follows that the constraint equations are The time derivatives of these functions have the expressions Then, the first order approximation of the the constraint equations in the neighborhood of the configuration (R 1 , u 1 , R 2 , u 2 ), has the expression , and

Spherical hinge
For a spherical hinge we start with the conditions or, with the notation used in the previous section, Then, the constraint function is with derivatives with respect to time or, in terms of angular velocities, Thus, the corresponding form of the C operator is 3.3.4Door hinge as a point slider and a spherical hinge The operators corresponding to a point slider and to a spherical hinge can be composed together to give that of a door hinge:

Integration of the kinematic-compatibility equations
Let E be the Euclidean group.The elements of E are the pair (R, u) with R ∈ Orth + the rotation and u ∈ V the translation.The group operation on E is defined by the inverse of (R, u) is (R −1 , −Ru), and the unit element of the group is the pair (I, 0 ).This make E the semidirect product of Orth + and V .
The Lie algebra e of E is given by the semi-direct product of the Lie algebras Skw and V .The exponential map exp : e → E transforms the pair (W , w ) ∈ e in the pair (R, u) ∈ E, with here θ denotes the norm of the axial ω of W and A the operator Let G be the direct product group of In components, the elements g of G are written g = (g 1 , g 2 , . . ., g N ), with g i = (R i , u i ) ∈ E. The i-th component g i of g defines the rigid transport of the i-th body.Then the Lie algebra g of G is the direct product of n copies of the Lie algebra e of E, and the exponential map exp : g → G is the map that transforms (u 1 , u 1 , . . ., u N ) ∈ g in g = (g 1 , g 2 , . . ., g N ) ∈ G, with g i the image of u i under the exponential map of the Euclidean group E.

Newton-Raphson method on R n
To introduce the notation adopted, we briefly recall here the well-known Newton-Raphson method on R n .Given a map f : R n → R m , we consider the problem of solving the nonlinear equation f (x) = 0 .(20) The iteration of the Newton-Raphson method for solving this equation is defined by Let p xn the function defined by p xn (v ) = x n + v , and let f the composed function as ∇p xn = I; in particular, for v = 0, we have It follows that the Newton-Raphson iteration can be rewritten in the form

Newton-Raphson Method on Lie groups
Let G be a Lie group.For every h ∈ G, the right translation is the map Given a map f : G → R m , we consider the problem of solving the nonlinear equation f (x) = 0 , (28) in which x is the collection of rotations and translations of all bodies, (R i , u i ), i = 1, . . ., N , and f is collection of all constraint functions c that we detailed in the previous subsections for each type of constraint.
Let f be the composed function because the exponential of 0 ∈ g is the identity e ∈ G, and the differential of the exponential map at zero is the identity of g, we have It follows the Newton-Raphson iteration on Lie group in which δv is the collection of all the increments (δω i , δw i ), i = 1, . . ., N .

Folding example
As the present work is mainly focused on the study of isostatic thick origami structures, we present just one folding example.We apply the above-described procedure to fold the triangular portion of Yoshimura pattern shown in Figure 2 (a).This triangulated origami structure is composed by isosceles triangles with a long side of 1.8 units and two equal angles of 22.5 degrees.The assembly has P = 30 panels with H = 35 hinge lines.The number of edges on the boundary is E b = 20.Hence, the number of internal independent mechanisms is M = E b − 3 = 17.We consider all hinges to be door hinges.The folding process is prescribed by imposing the relative angular velocities at the 15 vertical hinges to be equal to 0.005 rad/s, while maintaining vertex 6 fixed in space, the vertices 2, 4, 8, 10, 12, 13, 14, 15, 16 moving only along the x axis, and the motion of vertex 14 blocked along the x and y direction.
Figure 2 (b) shows the configuration reached after 200 s of simulation time.Figure 2 (c) shows snapshots of the folding process taken every 100 s of simulation time.

Static analyses
In this section, we first list the internal constraint reactions for each type of constraint, and then we report our analyses regarding the Miura-ori, the Yoshimura, and the Kresling crease patterns.

Internal constraint reactions
We recall that in Section 3, the first-order approximation of the internal constraint equation were obtained in the form c 0 + C δ = 0 , with c 0 the constraint functions, δ the increment of angular and translational velocities, and C an operator whose form depend on the considered constraint.
By assuming the constraints to be ideal, the vector r = (m 1 , f 1 , m 2 , f 2 ) of the reactive action at hinges, with m i , f i , i = 1, 2, the internal couple and the internal force exchanged at a hinge and reduced to the origin of the reference frame, is given by with λ the vector of internal actions expressed as Lagrange multipliers and B = C T .We detail below the form of the B operator for each constraint type.For a sliding hinge, from (4), we have λ = (λ 1 , λ 2 , λ 3 , λ 4 ) and for a sliding hinge as a double point slider, from (13), we have instead For a door hinge, from (6), we have λ = (λ 1 , λ 2 , λ 3 ) and Fig. 3 Two panels connected by a door hinge (a) and by a sliding hinge (b).Representations of the internal actions applied to a panel edge.The parallel moment M P is always null; the parallel shear force T P is null for sliding hinges.
for a door hinge as a point slider and a spherical hinge, from (15), we have instead The equilibrium equations of the whole origami structure are obtained by imposing, for each body, the resultant of all the internal constraint reaction and the external loads applied to it, and the resultant moment of the same with respect to a fixed pole, to be null.

Examples
In order to report the results of the static analyses, we reduce the constraints reactions to the mid points of the corresponding edges shared by adjacent panels, and project them along the directions of a right-handed local frame {n, t, m} attached to one of the two panels, with n the outward normal to the edge in the panel plane, t the tangent to the edge, and m the normal to the panel.In this way, the internal force and couple applied to a panel by the adjacent one are described by the triplets (N, T P , T O ) and (M T , M P , M O ), respectively, with N the normal force, T P the parallel shear force, T O the orthogonal shear force, M T the twisting moment, M P the parallel moment, and M O the orthogonal moment, as illustrated in Figure 3.The parallel moment is always null, while the parallel shear force is null for sliding hinges.
We introduce our gallery of examples by describing the simple case of a panel ring formed by four rectangular panels connected by four door hinges, which is shown in Figure 4 (a).The kinematics of this assembly is analogous to that of a planar fourbar linkage and therefore there is just one internal mechanisms, M = 1.Since there are no internal vertices, the number of self-stress states is S = 3, according to both the bar-framework model and the panel-hinge model.In addition, it is easy to check by straightforward equilibrium calculations that the internal actions N and T O must be null for all panel hinges.Then, by taking as parameters the parallel shear force, the twisting moment, and the orthogonal moment at one hinge, one can compute explicitly the internal actions for the three self-stress states, which are represented Fig. 4 A four-panel origami ring (a) and a representation of its three self-stress states (b,c,d).
in Figure 4 (b,c,d).It is easy to check that panel rings with more than four panels, P > 4, always have S = 3 and M = P − 3.
It is worth noticing that by replacing door hinges by sliding hinges it is only possible to eliminate the self-stress state with nonzero parallel shear T (Figure 4 (b)).The introduction of more than one sliding hinge would generate additional internal mechanisms without affecting the other two independent self-stress states.

Miura-ori patches
Similar to the case of the four-panel ring, the single-vertex Miura-ori origami shown in Figure 5 (left) has M = 1 and S P H = 3, when all the four hinges are door hinges.However, by performing numerical calculations with the proposed formulation, we checked that at non-singular configurations the replacement of any three door hinges with three sliding hinges makes the degree of static indeterminacy equal to zero.
Regular Miura-ori origami structures always possess just one internal mechanisms, no matter the number of panels they are composed of.Therefore, the degree of static indeterminacy always increases with the size of the Miura-ori patch.We found in our analyses that that there is a class of Miura-ori patches whose degree of static indeterminacy can be made equal to zero.The elements of this class can be constructed . . .Fig. 5 Left: a single-vertex Miura-ori origami with null degree of static indeterminacy.There are three sliding hinges and one door hinge (dashed line).Right: sequential assembly of isostatic Miuraori origami.There is just a single door hinge (dashed line), all the other hinges are sliding hinges.At each step, two panels and three sliding hinges are added.
sequentially, by starting from a single-vertex Miura-ori analogous to the one just described, and by iteratively adding two panels and three sliding hinges, as indicated in Figure 5 (right).We now give an example of calculation of internal actions for the case of the cantilevered Miura-ori beam shown in Figure 6.Each panel has the shape of a right trapezoid with height of 2 units, and bases of 1 and 2 units.The configuration in Figure 6 (b) is obtained by folding the initially flat configuration in Figure 6 (a) so that the angle between the normals to panels 2 and 3 is equal to 27.4 degrees.Table 1 reports the corresponding vertex coordinates.In the configuration of Figure 6 (b), two corner vertices on the short side of the assembly are pinned to the ground, and the middle vertex on the same side is constrained not to move along the longitudinal direction.At the other end of the assembly, a vertical point loads of unit magnitude is applied.The assembly has one door hinge (dashed line in Figure 6), and the remaining hinges are sliding hinges.The resulting internal and external constraint reactions are shown in Table 2 and Table 3, respectively.

Yoshimura wedges
We pass now to the Yoshimura origami structure described in Section 3.4.3.We consider the configuration shown in Figure 2 (b) and reassign vertex constraints as follows.Al base vertices, those on the plane z = 0, are fixed in the y and z direction, with vertex 6 constrained also in the x direction.We checked that with this set of vertex constraints, there are no internal mechanisms, M = 0.This origami structure has 6 internal vertices.When all hinges are door hinges, the degree of static indeterminacy is S P H = 3V i = 18.We verified that by replacing three door hinges at each internal vertex with three sliding hinges (cf. Figure 7) a isostatic structure is obtained.We considered the case of a vertical unit load applied to the tip of the structure, as shown in Figure 7. Table 4 lists the nodal coordinates of the analyzed configuration, while Table 5 and Table 6 reports the internal and external reactions, respectively.We checked that it is possible to obtain similar isostatic Yoshimura wedges with any number of base elements in an anologous way by constraining all base vertices in the transversal (y) and vertical (z) directions, with one of them constrained also in the folding direction (x), and by replacing three door hinges with sliding hinges at each internal vertex.

Kresling columns
We now pass to consider multimodular Kresling columns with triangular base.Each module of a Kresling column is composed by P = 6 triangular panels connected by H = 6 door hinges in a three-fold cyclic-symmetric configuration, as shown in Figure 8.The bottom and top horizontal bases of a module appear rotated relatively to each other by a certain twist angle θ ∈ (π/3, π).We remark that given the panels of a Kresling module, it is possible to assemble them into two distinct configurations with opposite handedness, mirror image of each other.
A one-module assembly has no internal vertices, and we found that is always isostatic (M = 0, S = 0), except for values of the twist angle equal to π 6 , or to − 5 6 π, at which the configuration is singular and admits one internal mechanism and one self-stress state (M = 1, S = 1).Figure 8 (a) shows the configuration of a Kresling module with height 20 units, bases inscribed in a circle with diameter of 20 units, and twist angle θ = π 6 .The corresponding vertex coordinates are listed in Table 7.We y x z Fig. 7 A Yoshimura origami wedge shaped as a cantilevered arch loaded at the tip.Door hinges are denoted by dashed white lines, the rest of them are sliding hinges.
computed the values of the internal actions in the self-stress state reported in Table 8. Figure 8 (b) shows a (isostatic) configuration with same height and base dimensions, and θ = π 9 , simply supported at the bottom base, and subjected to three vertical unit loads applied to the top vertices.Table 9 lists the vertex coordinates, while Table 10 and Table 11 reports the internal and external reactions, respectively.
A two-module column obtained by juxtaposition of two modules with same dimensions and opposite handedness is shown in Figure 9 (a).In addition to the hinges belonging to each individual module, three more hinges realize the junction between the two module.This origami structure possesses three internal vertices and therefore has S P H = 9 degrees of static indeterminacy in non-singular configurations, while M = 0. We checked that by realizing all hinges as sliding hinges, except for those belonging to one of the two modules, realized as door hinges, the assembly is isostatic.Isostatic multimodular Kresling columns can be obtained in a similar way. Figure 9 (b) shows the case of a three-module tower.Tables 12 and 13 shows the internal actions of the simply supported structures in Figure 9

Conclusions
We presented a design strategy aimed at obtaining isostatic thick origami structures by making use of sliding hinges to replace conventional hinges.This strategy can effectively reduce the degree of static indeterminacy, and in many cases it can successfully make certain types of origami structures isostatic.We provided exact and linearized kinematic constraint equations and a numerical procedure to simulate the folding process by controlling the relative angular velocities at some creases of the origami.Equilibrium equations relating the external loads to the internal actions exchanged by adjacent panels were obtained by duality.We analyzed a series of noteworthy examples of isostatic thick origami structures.We showed that Miura-ori crease patterns in a certain class admit an isostatic realization.We found that Yoshimura wedges can always be made isostatic.Finally, we demonstrated isostatic multimodular Kresling origami.The present method can be employed to discover, investigate, and design other types of isostatic thick origami structures and aid in the determination of self-stress states when this is not possible.

Fig. 1
Fig. 1 The square twist origami (a) and the corresponding bar framework idealization (b).There are P = 9 panels, H = 12 hinge edges (solid light gray lines), and E b = 8 boundary edges (solid black lines).Each of the five quadrilateral panels is triangulated with one additional edge (dashed light gray line) and transformed into a block by adding one more edge (dashed black line).There are V i = 4 internal vertices (black dots).The bar framework model of this origami has S BF = 1 selfstress state, and the panel-hinge model has S P H = 13, while M = 1 in both models.

Fig. 2
Fig. 2 Folding of a triangular portion of Yoshimura origami.(a) Initial configuration in the x − z plane.(b) Partially folded configuration.(c) Snapshot of the folding process.
(a) and (b), respectively, resulting from a vertical loading by three unit forces on the top vertices.

9 Fig. 8
Fig. 8 Two different versions of a Kresling module with triangular base.All hinges are door hinges in both (dashed white lines).Configuration (a) is singular with one infinitesimal mechanism and one self-stress state; configuration (b) is isostatic.Configuration (b) is simply supported on the ground, and subjected to vertical unit forces.

Fig. 9
Fig. 9 Two isostatic Kresling towers with two (a) and three (b) modules.Six hinges on one of the modules are door hinges (dashed white lines).

Table 8 Table 9
Internal actions in the self-stress state of the Kresling module in singular configuration shown in Figure 8 (a)).Only the values at two hinges are shown, the remaining hinges feature the same values according to the three-fold cyclic symmetry of the assembly.-02 -3.5127e-02 -3.8808e-02 -1.3086e-02 Vertex coordinates for the Kresling module in Figure 8 (b).

Table 1
Nodal coordinates for the Miura-ori cantilever in Figure 6(b).

Table 2
Internal actions in the origami structure in Figure 6 (b).

Table 3
Reaction forces at supports, and panels which they are applied to, for the origami structure in Figure6(b).

Table 5
Internal actions in the origami structure in Figure7.

Table 6
Reaction forces at supports, and panels which they are applied to, for the origami structure in Figure7.

Table 7
Vertex coordinates for the Kresling module in singular configuration in Figure8(a).

Table 10
Internal actions for the Kresling module in Figure 8 (b).Only the values at two hinges are shown, the remaining hinges feature the same values according to the three-fold cyclic symmetry of the assembly.

Table 11
Reaction forces at supports, and panels which they are applied to, for the Kresling module in Figure 8 (b).

Table 12
Internal actions for the Kresling two-module assembly in Figure9 (a).Only the values at five hinges are shown, the remaining hinges feature the same values according to the three-fold cyclic symmetry of the assembly.

Table 13
Internal actions for the Kresling three-module assembly in Figure 9 (b).Only the values at eight hinges are shown, the remaining hinges feature the same values according to the three-fold cyclic symmetry of the assembly.