Folding Simulation of Rigid Origami with Lagrange Multiplier Method
FFolding Simulation of Rigid Origami with LagrangeMultiplier Method
Yucai Hu a , Haiyi Liang a,b, ∗ a CAS Key Laboratory of Mechanical Behavior and Design of Materials, Department ofModern Mechanics, University of Science and Technology of China, Hefei, Anhui 230026,China. b IAT-Chungu Joint Laboratory for Additive Manufacturing, Anhui Chungu 3D printingInstitute of Intelligent Equipment and Industrial Technology, Wuhu, Anhui 241200, China.
Abstract
Origami crease patterns are folding paths that transform flat sheets into spatialobjects. Origami patterns with a single degree of freedom (DOF) have creasesthat fold simultaneously. More often, several substeps are required to sequen-tially fold origami of multiple DOFs, and at each substep some creases fold andthe rest remain fixed. In this study, we combine the loop closure constraintwith Lagrange multiplier method to account for the sequential folding of rigidorigami of multiple DOFs, by controlling the rotation of different sets of creasesduring successive substeps. This strategy is also applicable to model origami-inspired devices, where creases may be equipped with rotational springs andthe folding process involves elastic energy. Several examples are presented toverify the proposed algorithms in tracing the sequential folding process as wellas searching the equilibrium configurations of origami with rotational springs.
Keywords: rigid origami, Lagrange multiplier, loop closure constraint,sequential folding, elastic folding
1. Introduction
Origami is the art of folding a sheet of paper into artistic three-dimensionalobject. Origami artists have shown that interesting shape morphing can beattained by folding along crease pattern. Recently, origami goes beyond artsand inspires engineers and scientists to design novel mechanisms, structures, and ∗ Corresponding author
Email address: [email protected] (Haiyi Liang)
Preprint submitted to Journal of L A TEX Templates June 11, 2020 a r X i v : . [ c s . C E ] J un aterials such as deployable space solar sails [1, 2], medical stents [3], fold coresof sandwich structures [4, 5], foldable robots [6], meta-materials [7, 8, 9], etc.To facilitate understanding and design of origami and its inspired derivatives,parametric equations have been derived for regular origami tessellations madeup of identical unit cells [10, 9, 11, 12], which would become tedious and evenintractable for irregular crease patterns. The finite element method is applicablein general with the facets modeled by the shell elements and creases by elastichinges. It can provide detailed elastic deformation of the facets at the cost ofmodelling and computational time, which may not be of major concern. Efficientcomputational approaches are required to capture the global deformation of theorigami.The bar-and-hinge model can be regarded as a simplified finite elementmethod to account efficiently for folding kinematics when origami is subjectedto external forces or torques [13, 14, 15, 16, 17]. In this model, all crease linesare represented by elastic springs, which allow in-plane stretching as well as out-of-plane bending of facets. Schenk and Guest first introduced the bar-and-hingemodel for the mechanical analysis of origami where the origami structure is rep-resented by a pin-jointed truss framework, assuming infinitesimal deformation[13]. Fuchi et al. combined this bar-and-hinge model with topology optimiza-tion techniques for the design of origami-based mechanism [18]. Filipov et al.presented several improvements, including new triangulation schemes for thequadrilateral facets and stiffness parameters obtained from the sheet material,for realistic modeling of origami [14]. Recently, to solve large deformation prob-lems involving bifurcations and multistability, nonlinear bar-and-hinge modelshave been developed which synthesize techniques including nonlinear bar ele-ments, rotational springs with penalty near the fully folded state to avoid localpenetration, solvers of arch-length method, schemes for tracing the bifurcations,etc [15, 16]. Bar-and-hinge model can also be used to simulate the rigid origamifolding by projecting the motion into the nullspace of the global constraint ma-trix, in the limit of infinite large stiffness of bar elements [13, 19].For origami with stiff facets connected by soft creases, it would be muchharder to bend the facets than to fold along the creases. The rigid origamimodel is thus proposed in which the facets are assumed rigid and only rotationsare allowed along the creases. The condition of loop closure constraint has tobe fulfilled at every interior vertex for all the rigid facets to be compatible witheach other during the folding process [20, 21]. Wu and You have investigated thefolding of rigid origami based on the rotating vector model which describes the2oop closure constraint using quaternions [22]. For a given crease pattern, thethree-dimensional (3D) folded form of rigid origami can be uniquely determinedby fold angles which are the supplementary angles of the dihedral angles betweentwo adjacent facets, see Fig.1. The kinematic folding of the rigid origami canbe described by the parametric equations for origami tessellations made upof identical unit cells or origami with a single vertex. Wei et al. calculatedanalytically the Poisson’s ratio and stiffness of the Miura-ori tessellation basedon the geometry of the unit cell [10]. Assuming symmetry between the foldangles, Hanna et al. investigated the kinematics of degree-8 Waterbomb baseof single vertex with symmetric 8 creases [12]. Chen et al. presented thoroughkinematic folding analysis for degree-6 Waterbomb base of both thin and thickorigami [23].Though rigid origami model has been widely accepted in the literature dueto its simplicity, most of the previous works are elaborated on the origami tes-sellations consisting of identical unit cells or the origami base. To deal with rigidorigami of irregular crease pattern, Tachi developed a numerical method basedon a modified version of the loop closure constraint developed by belcastro andHull [20], to simulate the folding motion of rigid origami. In this algorithm, theinfinitesimal rotations of facets were calculated by projecting fold angle changesinto the linearized constraint space with the numerical residual compensatedby a single Newton-Raphson iteration [24, 25]. However, it may be seen thatthe folding increment and numerical residual are not strictly controlled for eachfolding step.Rigid origami model can be extended to search the elastic equilibrium config-urations under the competition between rotational springs mounted at creaseswith different rest angles, for example, in the case of self-folding mechanisms androbots with smart material actuators mounted at the folds [26, 27, 28]. Bruncket al. derived the covariant energy and the associated geometric constraint (i.e.,constant sector angles) for a vertex with n creases using the unit vectors alongthe creases as the variables [29]. Wang and Qiu adopted pseudo-folds to approx-imate the bent configuration in which the nodal coordinates and the pertinentconstraints are expressed with respect to the fold angles [30]. Both the afore-mentioned two studies focus on the rigid origami of single vertex. The groundstructure, a potential structure equipped with a sufficient number of pseudo-folds, has been employed for the origami design by considering the stiffness ofthe crease rotational springs as the weight function [31].Most of existing studies focus on origami of single DOF, such as Miura-ori3hat folds simultaneously. In many cases, crease pattern is of multiple DOFsand usually requires sequential folding, so that it takes several distinct steps tofold an origami and at each substage some creases fold and the rest remain fixed.This complicates the folding analysis, and leaves analytical solution out of reach.In this paper, we pay attention to sequential folding analysis of rigid origamiwhen the folding sequence is known. The proposed algorithm is presented forthe folding analysis of the rigid origami without inner holes. In the followingsections, the rigid origami model with the fold angles as the variables is reviewedfirst. Then, the Lagrange multiplier method is applied to control a subset ofcreases which drives the folding of the rigid origami. The numerical residual ofeach folding step is eliminated by the Newton-Raphson method, and thus thefold angles yield valid configuration. In addition, by considering the rotationalsprings at the creases, an algorithm based on the Lagrange multiplier methodis presented to search the equilibrium configuration of the rigid origami. It isshown that the projection method by Tachi [24] agrees with a special case of thecurrent physical model. The algorithms are then verified by several examples.
2. Rigid origami model
The rigid origami model of Tachi [24, 25] is reviewed in this section for com-pleteness. The origami is commonly designated by the crease pattern consistingof vertexes and creases. Vertexes are points on the origami paper and eachcrease is a line joining two neighboring vertexes along which origami paper isfolded. Fig.1(a) shows an isometric view of a typical vertex P with n creasesjoining it on an unfolded or flat origami paper. The creases are numbered from1 to n anticlockwise. Creases i and i + 1 define the sector angle θ i . It is trivialthat the sum of the n sector angles is 2 π . Before folding, AOB is a straight lineperpendicular to crease i on the origami paper, i.e., ∠ AOB is π as shown in thefigure. After folding, ∠ AOB would be denoted as π − ρ i in which ρ i is the foldangle of crease i . For the rigid facets to be compatible with each other, the loopclosure constraint around the vertex P is F ( ρ P ) = χ , χ , · · · χ n, = I . (1)4 a) (b) θ i- θ θ i- π ABO θ n n iθ i xyz n- θ n- i+ i- θ i- π-ρ i AB O i θ i i+ L i- M i- N i- L i M i N i P Figure 1: (a) A vertex with n creases before folding. θ + θ + · · · + θ n = 2 π . AOB is astraight line perpendicular to crease i with both A and B on the origami and ∠ AOB = π . (b)AOB is kinked after folding. ∠ AOB becomes π − ρ i , where ρ i is the fold angle of crease i .The triplet, L i , M i and N i , defines a local Cartesian coordinate system for the sector facet( i, i + 1): L i is the unit vector along crease i ; N i is the unit normal of the sector facet ( i, i + 1)and M i = N i × L i where “ × ” is the cross product. The indexes ( i -1) and ( i +1) are to beinterpreted cyclically. where χ i − ,i = cos θ i − − sin θ i − θ i − cos θ i −
00 0 1 ρ i − sin ρ i ρ i cos ρ i ; (2) I m is the m × m identity matrix and ρ P = { ρ , ρ , · · · , ρ n } is the vector con-taining the n fold angles. The derivation of Eq.(1) is given in Appendix A. Asthe facets bounded by adjacent creases are rigid, the sector angles are constantsand Eq.(1) are nonlinear constraints on the fold angles. It can be shown thatthe matrix ∂ F ∂ρ i is antisymmetric for compatible fold angles (see Appendix A forthe proof), i.e., ∂ F ∂ρ i (cid:12)(cid:12)(cid:12)(cid:12) ρ P = − c i b i c i − a i − b i a i for i = 1 , , · · · , n . (3)For each folding step in simulating the folding process, the iteration startswith an approximate solution ρ iP = { ρ i , ρ i , · · · , ρ in } and the associated residualis R ( ρ iP ) = F ( ρ iP ) − I . (4)5he aim of iteration is to find an increment or infinitesimal folding, ∆ ρ P , suchthat the residual for ρ i +1 P = ρ iP +∆ ρ P is eliminated. Substitute ρ i +1 P into Eq.(4)and expand the residual about ρ iP as, R ( ρ i +1 P ) = F ( ρ iP + ∆ ρ P ) − I (cid:39) F ( ρ iP ) + n (cid:88) j =1 ∂ F ∂ρ j (cid:12)(cid:12)(cid:12)(cid:12) ρ iP ∆ ρ j − I . (5)Thus, n (cid:88) j =1 ∂ F ∂ρ j (cid:12)(cid:12)(cid:12)(cid:12) ρ iP ∆ ρ j = −{ F ( ρ iP ) − I } = − R ( ρ iP ) . (6)It is assumed that the increment is infinitesimal and ρ iP is infinitesimally per-turbed from the valid state, thus the antisymmetry of the derivative matrix isapproximately retained. Thus, only three of the nine equations in Eq.(6) areindependent. Let a j , b j and c j be, respectively, the entries with the indexes(3,2), (1,3) and (2,1) in ∂ F ∂ρ j (cid:12)(cid:12)(cid:12) ρ iP , the independent equations of (6) can be writtenas a a · · · a n b b · · · b n c c · · · c n ∆ ρ ∆ ρ ...∆ ρ n = − F , F , F , (7)where F i,j is the ( i, j )-entry of the matrix F .For general crease pattern consisting of N V i interior vertexes, N Ei creases(interior edges) and no holes, Eq.(7) can be formed for every interior vertexwith the incident fold angles. Thus, there are 3 N V i constraints on the N Ei foldangles. Let ρ i be the vector of approximate fold angles renumbered globally,Eq.(7) for every interior vertex can be collected and expressed as C ∆ ρ = − r (8)where C is the global linearized constraint matrix of dimension 3 N V i × N Ei ; ∆ ρ is the vector of fold angle increment and r is the vector of residuals. Note thatboth C and r are evaluated at ρ i . If ρ i is a vector of valid fold angles, r is zeroand Eq.(8) reduces to the linearized constraint on the infinitesimal fold angles.Starting with a valid state and the vector of intended fold angle changes ∆ ρ ,Tachi introduced the Euler method to simulate the folding process by projecting∆ ρ into the nullspace of C and compensating the accumulated numerical error624]: ∆ ρ = [ I N Ei − C + C ]∆ ρ − C + r (9)where C + denotes the Moore-Penrose pseudoinverse of C . However, it is hardto control the folding increment using the above projection method.
3. Sequential folding using Lagrange multiplier method
Due to the loop closure constraint around every inner vertex, the foldingmotion of the origami can usually be driven by folding a subset of creases. Forinstance, by controlling anyone of the fold angles of the Miura-ori fold, the entireorigami can be folded simultaneously since the structure is a mechanism witha single DOF. For general origami, the DOF of the origami in the partiallyfolded state can be obtained from the nullity of the linearized constraint matrix C . Let I c be the set of controlled fold angles whose values are prescribed,i.e., ρ j = ¯ ρ j for j ∈ I c . The folding process consists of a number of foldingsteps. For the i -th folding step, the increment for the controlled angle ρ j is∆ ρ j = ¯ ρ i +1 j − ¯ ρ ij = f j where f j is self-defined. Except the controlled fold angles,the other components in ρ i are infinitesimally perturbed from the valid state.A direct approach is to eliminate the controlled fold angles from Eq.(8) which,however, alters the structure of the system matrix and can be inconvenientwhen different sets of creases are controlled in multiple folding stages. In thefollowing, the Lagrange multipliers λ = { λ j } are introduced for the controlledfold angles and the functional is Π (∆ ρ , λ ) = 12 ( C ∆ ρ + r ) T ( C ∆ ρ + r ) + (cid:88) j ∈ I c λ j (∆ ρ j − f j ) . (10)The solution is given by the point where the derivatives with respect to ∆ ρ and λ are zero: (cid:34) C T C AA T (cid:35) (cid:40) ∆ ρλ (cid:41) = (cid:40) − C T rf (cid:41) (11)in which A = [ e i ] i ∈ I c with e i being a column vector of length N Ei with 1 in the i -th position and 0 in every other position and f = { f j } . As the system can be7nder-determined, the minimum length solution for Eq.(11) is (cid:40) ∆ ρλ (cid:41) = (cid:34) C T C AA T (cid:35) + (cid:40) − C T rf (cid:41) . (12)It should be remarked that the vector of updated fold angles ρ i + ∆ ρ does notsatisfy the nonlinear loop closure constraint in Eq.(1) exactly in general. Byupdating C and r , and replacing f by so as to keep the controlled fold anglesunchanged, Eq.(12) is looped to reduce the residual r until that the convergencetolerance is satisfied for the i -th folding step. The procedure for a single foldingstep is summarized in Algorithm 1 which can be repeated to obtain the wholefolding motion.For each folding step, the resultant fold angles are compatible and the 3Dfolded form of the origami can be visualized by calculating the coordinates ofthe vertexes based on the crease pattern and fold angles. The procedure forcalculating the 3D folded form is presented in the Appendix B. Many origamiartworks involve a sequence of folding stages or sequential folding. With allthe creases and the folding sequence specified a priori , the simulation can beprocessed by updating I c and f for each stage. More details are exposed in theexample 5.2 by considering the sequential folding of a crane. Algorithm 1
A single folding step driven by controlled creases
Require:
Sector angles at all vertexes, θ ; Current fold angles, ρ i ; The specifiedincrement for ∆ ρ j with j ∈ I c , f ; Tolerance for the residual error, (cid:15) = 10 − . Ensure:
The vector of fold angles for the next step, ρ i +1 . Calculate C and r at ρ i ; Calculate the increment ∆ ρ from (11); Update ρ i +1 = ρ i + ∆ ρ and r ; while ( (cid:107) r (cid:107) / (3 N V i ) ≥ (cid:15) ) do Update C at ρ i +1 ; Solve ∆ ρ from (11) with f = ; ρ i +1 ← ρ i +1 + ∆ ρ ; update r at ρ i +1 ; end while4. Elastic folding using Lagrange multiplier method Rigid origami can be elastically folded and multistable by mounting rota-tional springs at crease lines. In fact, various actuators such as shape memoryalloy and shape memory polymer [26, 28] have been employed to actuate the8olding of origami inspired structures. Considering rotational springs at thecreases, the elastic energy of the origami structure is U ( ρ ) = 12 N Ei (cid:88) i =1 k i ( ρ i − ˜ ρ i ) (13)where k i and ˜ ρ i are, respectively, the rotational spring stiffness and rest angleof crease- i and the former is k i = kL i in which L i is the length of the creaseand the constant k is the stiffness per unit length. When the fold angle ρ isincreased by ∆ ρ , the energy increment is∆ U = U ( ρ + ∆ ρ ) − U ( ρ ) = 12 ∆ ρ T H ∆ ρ + d T ∆ ρ (14)where H = diag. { k , k , · · · , k N Ei } and d = [ k ( ρ − ˜ ρ ) , k ( ρ − ˜ ρ ) , · · · , k N Ei ( ρ N Ei − ˜ ρ N Ei )] T = ∂U∂ ρ (15)which is essentially a vector of internal moments along the creases. The aim isto find ∆ ρ such that ∆ U is minimized while subjected to the loop closure con-straint of the rigid origami, see Eq.(1). By considering the linearized constrainton the increment in Eq.(8), the minimization of ∆ U subjected to nonlinear con-straints is transformed into a classical quadratic problem [32]. Introducing theloop closure constraints on internal vertices by the Lagrange multipliers λ , thefunctional on the infinitesimal increment can be written as∆ Π U = 12 ∆ ρ T H ∆ ρ + d T ∆ ρ + λ T ( C ∆ ρ + r ) . (16)Variations of (16) with respect to ∆ ρ and λ yields (cid:34) H C T C 0 (cid:35) (cid:40) ∆ ρλ (cid:41) = − (cid:40) dr (cid:41) . (17)The minimum length solution for Eq.(17) is (cid:40) ∆ ρλ (cid:41) = − (cid:34) H C T C 0 (cid:35) + (cid:40) dr (cid:41) . (18)In the case of rigid origami with triangular facets with no holes, it can be shownthat 3 N V i ≤ N Ei [33]. Assuming that all the linearized constraints in Eq.(8) areindependent, i.e., Rank( C ) = 3 N V i , explicit expression for the inverse matrix in9q.(18) can be derived and the increment is (see Section 16.2 of the book [32])∆ ρ = − ( H − − GCH − ) d − Gr with G=H − C T ( CH − C T ) − . (19)In the special case that all springs are of the same stiffness, i.e., k i = k for i = 1 , , · · · N Ei , the diagonal matrix H reduces to the identity matrix I N Ei multiplied by k . Thus, we have G = C T ( CC T ) − = C + and Eq.(19) can besimplified as ∆ ρ = − k ( I N Ei − C + C ) d − C + r . (20)As d is the gradient of the energy, Eq.(20) means that the vector of incrementis given by the fastest decreasing direction of the system energy projected tothe linearized constraint space. It can be seen that Eq.(20) agrees with theprojection method in Eq.(9) (also see Eq.(14) of the reference [24]) when theinternal moment d is treated as the vector of intended fold angle changes. Here,the physical meaning is more clear and Eq.(18) should be used for general cases.It should be remarked that ∆ ρ suggests the optimal direction to decreasethe system energy while satisfying the linearized closure constraint. The steplength of the searching should, however, be restricted such that the changesin the fold angles are still infinitesimal. In this light, a step length factor c isintroduced and ρ i +1 = ρ i + c ∆ ρ max(abs(∆ ρ )) (21)where max(abs(∆ ρ )) extracts the largest magnitude of the components in thefold angle increment vector. The factor enforce that the largest change in a foldangle is c during a single folding step and we restrict that c ≤ π/
36. Similar tothe discussion following Eq.(12) in Subsection 3, the fold angle ρ i +1 does notexactly fulfill the nonlinear loop consistency constraint in Eq.(1). The residual r is to be compensated with iterations of ∆ ρ = − C + r . Since the length ofthe fold step is kept small, the energy will typically decrease monotonously inthe initial folding steps. A fold angle which increases/decreases monotonouslybefore reaching the local energy minimum is chosen as the characteristic foldangle and indicated by ρ a . When the increment of the characteristic angleis reversed, i.e., ∆ ρ a ( ρ ia − ρ i − a ) <
0, the local minimum state is within theregion bounded by ρ ia and ρ ia + c ∆ ρ a / max(abs(∆ ρ )). To converge to the local10inimum, the step length factor c is divided by two, i.e., c ← c/
2, and thealgorithm is summarized in Algorithm 2.
Algorithm 2
Folding driven by the rotational springs
Require:
Sector angles at all vertexes, θ ; Initial fold angles, ρ ; Spring stiffnessfor the crease- i , k i ; Characteristic fold angle, ρ a ; Initial step length factor, c ; Tolerances, (cid:15) and (cid:15) . Ensure:
Folding states ρ i and the converged fold angles for the local minimum. i ← , c ← c , ρ ← ρ ; while c > (cid:15) and i < specified maximum increment number do i ← i + 1; Calculate C and r at ρ i ; Calculate the increment ∆ ρ from (18); if i > ρ a ( ρ ia − ρ i − a ) < then c ← c/ end if ρ i +1 ← ρ i + c ∆ ρ / max(abs(∆ ρ )) and update r ; while ( (cid:107) r (cid:107) / (3 N V i ) ≥ (cid:15) ) do Update C at ρ i +1 ; ∆ ρ ← − C + r ; ρ i +1 ← ρ i +1 + ∆ ρ ; Update r at ρ i +1 ; end while end while5. Simulations of sequential/elastic folding In this section, the algorithms proposed in the previous sections 3 and 4 arevalidated by four examples of controlled origami folding, where the advantagesof combining the loop closure constrains and Lagrange multiplier method aredemonstrated. The first two examples are dedicated to the sequential foldingalgorithm, i.e., Algorithm 1. In subsection 5.1, as an example of single DOF,the simultaneous folding of Miura-ori is achieved by solely controlling one creasein one simulation step, and the folding kinetics is compared with the analyticsolutions. Subsection 5.2 illustrates the sequential folding of origami crane, antypical origami example of multiple DOFs, using three folding substeps, wheredifferent sets of creases are successively controlled with Lagrange multiplier.The last two examples in subsections 5.3 and 5.4, folding simulations of Water-bomb and Waterbomb tessellation, show the elastic bistability or equilibriumconfiguration can be reached iteratively using Algorithm 2, when rigid origamiare equipped with rotational springs at creases.11 a) (b)(c)
W L H αab ρ ρ ρ ρ α ρρ ρρ
12 34 ρ ρ (d) Figure 2: (a) Parameters for the Miura-ori unit cell; mountain and valley creases are indicatedby the black solid and black dashed lines, respectively; (b) partially folded configuration ofthe unit cell, the fold angle of valley crease is positive, i.e, ρ >
0, whilst those of mountaincreases are negative; (c) crease pattern for the Miura-ori fold consisting of 3 × a = b = 1 and α = π/
3; (d) the dimensions of the folded form. In (c) and (d),the fold angle of the red crease is controlled with the Lagrange multiplier.
The Miura-ori unit is made up of four identical parallelograms characterizedby the parameters a , b and α , see Fig.2(a) and (b). The folding motion of aMiura-ori fold with 3 × ρ asshown in Fig.2(c). For the planar state with all fold angles equal 0, the thirdequation in Eq.(7) will degenerate as c i , the (2,1)-entry of the derivative ma-trix in Eq.(3), equals 0 [25]. This reflects the fact that several folded formsare permissible as the mountain and valley crease assignment is missing. Inthe simulation, the initial fold angles are prescribed with small values, for in-stance ± ◦ , whilst the positive and negative signs are assigned for the valleyand mountain fold angles, respectively. The numerical error can be eliminatedby the iterations in the Algorithm 1.During the simulation, the controlled angle ρ , see the red crease in Fig.2(c)and (d), decreased from 0 ◦ to − ◦ by − ◦ in every folding step. Fig.3(a)shows the ρ versus ρ which is in agreement with the analytical solution (the12 a) (b) - ρ (ᵒ)
30 60 90 120 150 18001234567 -3-2.5-2-1.5-1-0.50 L Poisson's ratio W - ρ ( ᵒ ) - ρ (ᵒ) Figure 3: (a) The fold angle ρ versus ρ ; the insets are frames for ρ = − ◦ , − ◦ and − ◦ ;(b) the length L , width W and Poisson’s ratio versus ρ . For each variable, the markers arethe results obtained form the simulation whilst the dash lines are the analytical solutions[9]. dash line): tan( ρ /
2) = cos( α ) tan( ρ /
2) [34]. The insets in Fig.3(a) are thefolded form with ρ = − ◦ , − ◦ and − ◦ . The dimensions of the foldedform, see Fig.2(d), are measured at the end of each folding step. Fig.3(b)shows the change of width and length of the 3 × ν LW = − (cid:0) dLL (cid:1) / (cid:0) dWW (cid:1) [10],is calculated numerically as ν hLW = − (cid:0) L i +1 − L i L i (cid:1) / (cid:0) W i +1 − W i W i (cid:1) where the subscript i and i + 1 correspond to the number of folding step. It can be seen that thenumerical predictions for Poisson’s ratio agree with the analytical solution andthe accuracy can be improved if a smaller folding step is used. This examplevalidates the effectiveness of Algorithm 1 which is of potential to trace exactlythe geometric change of general rigid origami. It is worth mentioning that thefold angle ρ is folded exactly to the specified angles during the folding processwhich presents challenge for the projection method in the reference [24]. When folding the classical origami crane from a blank sheet by hand, asequence of folding steps are needed to attain the final shape due to its multipleDOFs of the crease pattern. Some of the creases become blocked and inactive,i.e., ρ i = − π, π , as the folding proceeds. In the simulation, all the creasesare specified a priori as shown in Fig.4(a). Fig.4(b) shows the paper model forthe finial folded crane. It is worth mentioning that, unlike Randlett’s flapping13 a) (b) Beak
Tail Beak
Tail
Figure 4: (a) Creases and boundary edges for the crane; (b) Photo of the folded paper crane.The beak and tail of the crane are indicated by the circle and square, respectively. bird [17], the current crane is rigid-foldable.Through folding by hand, three successive folding stages are identified. Thecrease pattern for the first stage is shown in Fig.5(a) where the highlightedcreases are controlled. The gray creases are inactive in this stage, i.e., thepertinent fold angles equal 0 throughout this stage. Thus, except the centralvertex, other vertexes are of degree-4 and are also flat-foldable. For flat-foldabledegree-4 vertex, it is known that the fold angles of opposite creases are equal inmagnitude [34]. Besides, the symmetry with respect to the diagonal creases isassumed. Thus, the fold angles at all the controlled creases are equal and theyincrease from 0 to π by a small amount per folding increment, for instance 5 ◦ in the simulation. The controlled fold angles ρ is chosen as a representativeand its history is shown in Fig.5(d). The fold angles of the uncontrolled creasesare calculated by the Algorithm 1. At the end of the first folding stage whenthe controlled fold angles reach π , the paper is folded flat. Three frames for thefirst stage are shown as insets in Fig.5(d).The crease pattern for the second folding stage is shown in Fig.5(b). Apartfrom the controlled creases inherited from the first stage (the highlighted graycreases), the highlighted solid creases numbered from 2 to 6 are also controlledwhich drive the folding of the second stage. The fold angles at creases 2, 3 and4 decrease from 0 to − π whilst those at crease 5 and 6 increase from − π to 0,see Fig.5(d) for the history of ρ . Three frames of the folded form are shown bythe insets in Fig.5(d) among which the last one shows that the paper is againfolded flat. As shown in Fig.5(c), despite the controlled fold angles inherited14rom the first and second stages, four creases numbered from 7 to 10 are addedto be controlled which drive the folding of the last stage. Fig.5(d) shows twoframes for the last stage and the history of ρ . In the last two stages, it can beseen that more creases become inactive, i.e., more fold angles keep constant at − π, π . To transfer from one stage to the next, the fold angles should beexactly controlled to reach the planar state which is enforced by the Lagrangemultipliers on the controlled fold angles. Fig.6(a) shows the circular waterbomb base of unit radius consisting of eightalternative mountain and valley creases around a single vertex. It has beensuggested as a test bed for actuated origami systems since the waterbomb basehas multiple DOFs and exhibits bistable behavior [12]. For simplicity, the springstiffness per unit length is taken as 1 in the following.First, we consider the symmetric case in which all mountain folds are ofthe same rest angle and so are the valley folds. Thus, the equilibrated foldedforms should also be symmetric. Referring to Fig.6(b), (c) and using sphericaltrigonometry, the analytical solutions for the mountain and valley folds are[ ρ m , ρ v ] = (cid:34) θ − π, (cid:16) √ θ − − √ θ (cid:17) − π (cid:35) ≤ θ < π (cid:34) π − θ, (cid:16) √ θ − √ θ (cid:17) − π (cid:35) π ≤ θ ≤ π . (22)where θ is the angle between −→ OS and −→ OA . The analytical solutions, [˜ ρ m , ˜ ρ v ] =[ − π/ , . · · · ] obtained by substituting θ = 5 π/ − π, π/
2] to search for the stable state similar to thatshown in Fig.6(b) and [ − π/ , π ] for that shown in Fig.6(c); the two initial statescorrespond to the downward compactly folded state with θ = 0 and the upwardcompactly folded state with θ = 3 π/
4, respectively. The searching paths of theAlgorithm 2 are indicated by the “*”s and “ ◦ ”s in Fig.7. The folded forms of thebistable states are similar to those in Fig.6(b) and (c) and thus are not plotted.The predictions for the fold angles and energy of the equilibrium states are inagreement with the analytical solution.15 ρ ρ (a) (b)(d) (c) ππ /2 -π /2 -π folding increment Mountain Crease; Valley Crease; Inac!ve Crease;
Figure 5: Sequential folding of origami crane. (a), (b) and (c) show the controlled creasepatterns for the first, second and third folding substeps, respectively; The foldable creasesin solid and dash line are mountains and valleys respectively; the creases under control arehighlighted in orange color, and the light gray solid lines are the inactive creases. (d) plotsthe history of ρ , ρ and ρ which correspond to the fold angles at the labeled creases in (a),(b) and (c); the insets show the 3D folded forms during the folding substeps. a) )c()b( Figure 6: (a) The crease pattern for the waterbomb base. (b) and (c) are the folded form ofthe bistable states for the case with symmetric rest angles, i.e., all mountain folds are of thesame rest angle and so are the valley folds; θ is the angle between OS and OA. E ne r g y Downward statesUpward statesAnalytical θ ≈ 0.9010 (51.62 (cid:101) ) Energy ≈ 3.5820 θ ≈ 1.9635 (112.5 (cid:101) ) Energy = 0
Figure 7: The energy versus θ in the searching process for the bistable states of the casewith symmetric rest angles by the Algorithm 2. The “*” at θ = 0 (“ ◦ ” at θ = 3 π/
4) is forthe downward (upward) compactly folded state and used as the initial state for the searchingprocess. a) (b) (c) Initial state ( ρ π ) - -5 /6 -2 /3 - /2 - /3 - /6 E ne r g y Downward statesUpward states ρ ≈ -1.3981 Energy ≈ 0.2878
Initial state ( ρ = - π /2) ρ ≈ -0.7730 Energy ≈ 0.7222
Figure 8: The waterbomb base with the unsymmetric rest angles: (a) the searching pathsby Algorithm 2 for the downward and upward stable states are indicated by “*”s and “ ◦ ”s,respectively; (b) and (c) are respectively the folded forms for the downward and upward stablestates where the facet 0-1-2 with the red edges are kept stationary. Next, we chose the rest angles as ˜ ρ i = ( − i πi +1 for crease- i with i =1 , , · · · ,
8. Analytical solutions are not available for this case. Similar to thesymmetirc case, the two initial states, i.e., downward and upward compactlyfolded states, are used to search the bistable states. As the states in the search-ing process are not symmetric, the angle θ is not well defined. Instead, thesearching path of the Algorithm 2 is indicated by the energy against the charac-teristic fold angle ρ shown in Fig.8(a). The folded forms of the bistable statesare shown in Fig.8(b) and (c), respectively. This example considers the waterbomb tessellation made up of degree-6bases, see the square with red edges in Fig.9. The rotational stiffness per unitlength is taken as k = 1. The rest angles for the mountain and valley creases are,respectively, set as ˜ ρ v = − ˜ ρ and ˜ ρ m = ˜ ρ and three cases with ˜ ρ = π/ , π/ π/ ρ in-dicated in Fig.9 is used to characterize the folded states during the searchingprocess. Fig.10(a) shows the searching paths of Algorithm 2 for the three casesand the corresponding equilibrium configurations are in Fig.10(b), (c) and (d).It can be seen that the waterbomb tessellation curves into a tube for ˜ ρ = π/ ρ increases. The computational costs of the three cases are,18 a b ρ Figure 9: The crease pattern for the Waterbomb tessellation made up of 5 × a = b = 1. respectively, 1.8s, 2.3s and 2.5s for the Algorithm 2. As the residual of eachstep is negligible, the intermediate states of the Algorithm 2 are valid and canbe visualized; the snapshots for the states with ρ = 0 ◦ , . ◦ , . ◦ , . ◦ and157 . ◦ are shown in Fig.11(a) through (e), respectively.
6. Closure
In this paper, we propose algorithms for rigid origami folding analysis usingfold angles as variables. By combining the loop closure constraint with La-grange multiplier method, we are able to model sequential origami folding whenthe crease pattern is of multiple DOFs. The introduction of Lagrange multi-plier method allows us to control the fold angles of different sets of creases, sothat some creases fold and the rest remain fixed at successive substeps. Newton-Raphson method is adopted in the algorithm to eliminate the numerical residual,and the geometric features of the origami can be accurately traced during thesequential folding simulations. This strategy is also extended to model rigidorigami with rotational springs mounted at the creases, which involves elasticenergy cost and competition between creases during the folding process. To findthe equilibrium configurations of origami with elastic rotational springs of differ-ent rest angles, we construct a functional to minimize the elastic spring energywhile enforcing the loop closure constraint with Lagrange multiplier method.The two algorithms are applicable to general origami structures without For reference, the MATLAB implementation was run on a desktop with Intel(R)Core(TM) i7-6700 (8 cores, 3.41Gz). a) y -0.50.51 x z (b)(c) y x z (d) x y
02 0 z E ne r g y ρ ≈ 1.58 (90.7 (cid:101) ) Energy ≈ 11.23 ρ ≈ 2.37 (135.6 (cid:101) ) Energy ≈ 4.14 ρ ≈ 2.75 (157.5 (cid:101) ) Energy ≈ 1.10
Figure 10: The searching paths (a) and the equilibrium configurations (b), (c) and (d) for thecases with ˜ ρ = π/ , π/ π/ a) (b)(c) x -101 y z x y z -24 050 401 x y z -24 -35 0 010 22 xy -1 1 z (d) x y
02 0 z (e) Figure 11: Snapshots of the intermediate states for the case with ˜ ρ = 7 π/ ρ = 0 ◦ , . ◦ , . ◦ , . ◦ and157 . ◦ where ρ is indicated by the red line. Appendix A. Loop closure constraint
This appendix first derives the loop closure constraint; then the derivative ofthe constraint matrix with respect to fold angles is proved to be antisymmetricfor compatible fold angles. In Fig. 1(b), the triplet, L i , M i and N i , defines alocal Cartesian coordinate system for the sector facet i -( i +1) with L i being theunit vector along crease i , N i being the unit normal of the sector facet i -( i + 1)and M i = N i × L i where “ × ” is the cross product. It is clear that[ L i , M i , N i ] = [ L i − , M i − , N i − ] χ i − ,i (A.1)where χ i − ,i given in Eq.(2) for i = 1 , , · · · , n and ( i -1) is to be interpretedcyclically. Starting from the sector facet 1-2 and looping around the vertexanticlockwise, recursive usage of the transformation in Eq.(A.1) yields[ L , M , N ] = [ L , M , N ] χ , χ , · · · χ n, . (A.2)Thus, the loop closure constraint is obtained as F ( ρ P ) = χ , χ , · · · χ n, = I . (A.3)Since χ i,j ’s are orthogonal, FF T always yield the identity matrix and itsderivative with respect to any fold angle vanish. With Eq.(A.3), it is clear that2224] ∂ ( FF T ) ∂ρ i (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ P = (cid:34) ∂ F ∂ρ i F T + F ∂ ( F T ) ∂ρ i (cid:35)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) ρ P = ∂ F ∂ρ i (cid:12)(cid:12)(cid:12)(cid:12) ρ P + ( ∂ F ∂ρ i ) T (cid:12)(cid:12)(cid:12)(cid:12) ρ P = . In other words, the matrix ∂ F ∂ρ i is antisymmetric when the fold angles are com-patible. Besides, from Eq.(A.3) and (2), we have ∂ F ∂ρ i = χ , · · · χ i − ,i − ∂ χ i − ,i ∂ρ i χ i,i +1 · · · χ n, (A.4)where ∂ χ i − ,i ∂ρ i = cos θ i − − sin θ i − θ i − cos θ i −
00 0 1 − sin ρ i − cos ρ i ρ i − sin ρ i . Thus, the entries of ∂ F ∂ρ i , i.e., a i , b i and c i in Eq.(3), can be obtained analyticallyfrom Eq.(A.4). Appendix B. Folded form for given crease pattern and fold angles