Origami provides both inspiration and potential solutions to the fabrication, assembly, and functionality of various structures and devices. Kinematic modeling of origami-based objects is essential to their analysis and design. Models for rigid origami, in which all planar faces of the sheet are rigid and folds are limited to straight creases having only zeroth-order geometric continuity, are available in the literature. Many of these models include constraints on the fold angles to ensure that any initially closed strip of faces is not torn during folding. However, these previous models are not intended for structures with non-negligible fold thickness or with maximum curvature at the folds restricted by material or structural limitations. Thus, for general structures, creased folds of merely zeroth-order geometric continuity are not appropriate idealizations of structural response, and a new approach is needed. In this work, a novel model analogous to those for rigid origami with creased folds is presented for sheets having realistic folds of nonzero surface area and exhibiting higher-order geometric continuity, here termed smooth folds. The geometry of smooth folds and constraints on their associated shape variables are presented. A numerical implementation of the model allowing for kinematic simulation of sheets having arbitrary fold patterns is also described. Simulation results are provided showing the capability of the model to capture realistic kinematic response of origami sheets with diverse fold patterns.

Introduction

Until recently, the term origami has been associated primarily with the ancient art of folding paper [1]. In origami, a goal shape is achieved from an initially planar sheet exclusively through folding. In this context, a fold is any deformation of the sheet in which the in-surface distance between any two points in the sheet is constant and the sheet does not self-intersect [2,3].

Engineering advantages of origami-inspired structures and devices include compact deployment/storage capability [4], a reduction in manufacturing complexity [5,6], and the potential for reconfigurability [7,8]. Existing and potential applications of origami solutions to device and structural design problems include deployable structures for space exploration [912], electronic components with improved properties [1315], robotic components [16,17], foldable wings [18], cellular materials [19], metamaterials [2023], and shelters [24,25], among others [2628].

Rigid origami, the special case of origami for which the planar faces of the sheet are inflexible [29,30] has been studied in the past and remains an active subject [31]. Rigid origami has been utilized for the design of deployable structures and architectural constructions [30,3234]. Theoretical modeling and simulation of rigid origami structures permit understanding of their kinematic behavior and the development of computational tools for their design. Rigid origami has been modeled using diverse approaches [29,35]. For example, Belcastro and Hull [36,37] presented a model for rigid origami derived by representing the deformation associated with folding using affine transformations. Their model provides constraints on the fold angles allowing for valid rigid origami configurations as well as mappings between unfolded and folded configurations. Tachi developed the Rigid Origami Simulator [29,38] for the simulation of rigid origami that also considered a set of constraints on the fold angles analogous to those presented in Refs. [36,37]. Using a similar approach, Tachi also developed Freeform Origami [39] for the simulation and design of freeform rigid origami structures represented as triangulated meshes [40].

Alternatively, truss representations [41] have been used wherein the polygonal faces of the sheet are triangulated, each fold or boundary edge end-point is represented by a truss joint, and each fold and boundary edge is represented by a truss member. The configurations for which the displacements of the truss joints do not cause elongations of the truss members represent valid rigid origami configurations. Additional constraints that allow the triangulated polygonal faces to remain planar are also considered for these models.

The majority of origami modeling approaches and design tools to date are based on the assumption of creased folds (see Fig. 1(a) for an example) that are straight line segments in the sheet that, upon folding deformation, the sheet has zeroth-order geometric continuity (G0) at such lines (i.e., the sheet tangent plane may be discontinuous at these folds). Curved creased folds are also allowed in origami (see Refs. [4244]); nevertheless, the focus of this work is on rigid origami-based structures for which curved folds are not allowed because their folding deformation induces bending of the faces joined to such folds [43].

The idealization of physically folded structures as sheets having creased folds has been useful in the modeling and design of several origami-inspired structures in the past [20,23,32,45,46]. However, this simplification may not be appropriate for structures having non-negligible fold thickness or constructed from materials that do not provide sufficient strain magnitudes to generate the high curvatures required for a creased idealization. For these structures, the obtained folded regions may not be accurately represented as creases but rather as bent sheet regions having higher-order geometric continuity. These folded regions are referred to in this work as smooth folds (an example of a sheet with smooth folds is shown in Fig. 1(b)).

There have been past efforts to model and analyze surfaces that exhibit bent and creased folding. Origami-based bent and creased surfaces have been simulated using collections of developable surface subdomains [4751]. Such advancements allow for the realistic animation and rendering of curved folds and combinations of creases and bent regions. However, none of the aforementioned works [4751] have considered constraints on the geometry and deformation of the bent folded regions that are required to preserve rigid faces as in analogous rigid origami models, which are essential when fold intersections or holes are present in the sheet. In view of this, a novel model for the kinematic response of rigid origami-based structures having smooth folds is presented in this work. The proposed model for origami with smooth folds includes rigid origami with creased folds as a special case, hence being more general.

Modeling of origami-based morphing of plate structures having significant thickness at the fold regions requires the arbitrary order of continuity of smooth folds. The present model is also useful in the kinematic analysis of sheets folded via active material actuation, where the achievable curvature at the folds is limited by the maximum strain magnitude provided by such active materials [27]. Examples of these active material-based origami structures include liquid-crystal elastomer [52,53], shape memory alloy [5457], dielectric elastomer [58,59], and optically responsive polymeric self-folding sheets [6062], among others [27,63]. The different assumptions regarding strain distributions at the fold regions associated with the implementation of these various materials require the arbitrary order of continuity considered in this work.

The basic origami concepts and a review of the model for rigid origami with creased folds extended in this work are presented in Sec. 2. Section 3 presents the newly proposed model for origami with smooth folds. It includes the geometric description of smooth folds and the fold pattern, constraints required for valid configurations, and the numerical implementation of the model allowing for kinematic simulation of sheets having arbitrary patterns of smooth folds. Section 4 presents the simulation results of sheets with diverse fold patterns that demonstrate the model capabilities and the resulting realistic kinematic structural response captured by the model. Finally, Sec. 5 provides a summarizing discussion and concludes the paper.

Concepts and Review of Origami With Creased Folds

This section presents various concepts of origami and briefly reviews a model of rigid origami with creased folds for the purposes of developing the novel extensions proposed in this work. The notation used in this work is also introduced.

The origami modeling approach adopted and extended herein is largely based on Refs. [36,37]. The considered continuum body denoted as the sheet, and the shape variables associated with the creased folds in the sheet (e.g., fold angles) are first described. The layout of the folds in the sheet (i.e., the fold pattern) is then determined by vertices (start-points and end-points of the line segments coincident with the creased folds in a planar reference configuration) and their connectivity. After the geometric parameters of the fold pattern are determined, constraints on the fold shape variables required for valid rigid origami configurations are described. The set of configurations that satisfy the stated constraints comprises the constrained configuration space. Continuous motion of the sheet is achieved by continuously altering the fold shape variables such that any attained configuration is contained in the constrained configuration space (i.e., the motion of the sheet is a continuous path in the constrained configuration space).

In the present approach, the required fold pattern data (e.g., vertex coordinates and fold connectivity) are provided as inputs. Information on methods to design fold patterns for specific applications is provided in Refs. [3,46,6471] and the references therein. The simulation of the continuous motion of the sheet is executed by incrementally updating the values of the fold shape variables using input guess increments and subsequently applies corrections such that the resulting set of fold shape variables satisfies the proposed constraints. Here, input guess increments for the fold shape variables are arbitrary and not obtained via any specific method. Motion planning in origami for various applications is being addressed by various researchers, and the reader is referred to Refs. [7276] for such works.

To begin, the various concepts of origami previously introduced must be formalized. The considered continuum body is denoted as the sheet, which in this work is a three-dimensional, orientable, path-connected surface with boundary. A detailed description of a sheet for origami with creased folds is provided in Ref. [3]. For origami with creased folds, the sheet is divided into various faces that are connected by straight edges corresponding to the creased folds. Every face comprising the sheet has the same aforementioned characteristics of the sheet (i.e., they are three-dimensional, orientable, path-connected surfaces with boundary).

The orthonormal vectors ei3, i = 1, 2, 3, with e3:=e1×e2 form the basis {e1,e2,e3} that defines the fixed global coordinate system. The reference configuration of the sheet, denoted S0, is defined such that it is contained in the plane spanned by e1 and e2 with its faces not overlapping each other, except at their shared boundary edges. Although several applications that utilize origami concepts do not start from a planar, nonoverlapping configuration (e.g., Refs. [17,23,77]), the reference configuration S0 has such characteristics in this work to agree with conventional origami modeling approaches [36,37] and for the sake of simplicity. The configuration of the faces comprising S0 is denoted P0iS0, i=1,...,NP, where NP is the number of faces in the sheet (i.e., S0=i=1NPP0i). The side of S0 with normal e3 is selected as the positive side of the sheet. An example of a sheet in its reference configuration S0 is provided in Fig. 2(a).

A current configuration of the sheet is denoted St where the parameter t indicates the history of deformation from the reference configuration (t = 0) to a current configuration (t > 0). Examples of a sheet in the reference configuration and a current configuration are presented in Figs. 2(a) and 2(b), respectively. The configuration of the faces comprising St is denoted PtiSt, i=1,...,NP, i.e., St=i=1NPPti. As previously stated, the focus of this work is on rigid origami, which imposes the following constraints for the configurations St [36,37]:

Definition 2.1. Valid configuration:A valid current configurationSthas the following characteristics: (i) the faces undergo only rigid deformations (i.e., they are neither stretched nor bent), (ii) the sheet is not torn (initially joined faces remain joined), and (iii) the sheet does not self-intersect.

The only nonrigid body deformations of the sheet are thus achieved by rotating adjacent faces relative to one another along their connecting creased fold in such a manner that the sheet only attains valid configurations during such deformations. Therefore, the sheet is G0 continuous at the creased folds.

Any configuration of the sheet is fully described by the set of the shape variables associated with the creased folds. The only shape variable associated with a creased fold describes the relative rotation between the two faces joined by such a fold and is denoted as fold angle.

Definition 2.2.Fold angle:The fold angleθ̂i(t)is defined as π minus the dihedral angle between the positive sides of the two faces joined by the ith creased fold.

Schematics showing the concept of fold angle are provided in Fig. 3. It is noted in the previous definition that each fold angle is a function of the parameter t. Specifically, each fold angle θ̂i(t),i=1,...,NF (where NF is the number of creased folds in the sheet) is a continuous function with respect to t since the motion of the sheet must be continuous (see Chap. 11 of Ref. [3]). For the remainder of this section, the dependence of θ̂i on t is left implicit to simplify the notation.

To avoid self-intersection of any pair of faces joined by a creased fold, the value of the associated fold angle must be contained in the interval [π,π]. The vector θ̂NF is constructed by collecting the fold angles as follows: 
θ̂=[θ̂1θ̂NF]
(1)

Fold Pattern.

The layout of the creased folds in S0 is presented in a fold pattern, often referred to as crease pattern in this context of origami with creased folds [3,36,37]. To describe the fold pattern, the vertices are introduced:

Definition 2.3.Vertices:The vertices are the start-points and end-points of the line segments coincident with the creased folds inS0. Each vertex has an associated position vector denotedvjspan(e1,e2).

The number of vertices located at the interior of S0 is denoted NI, and the number of vertices located at S0 (the boundary of S0, see Fig. 2(a)) is denoted NB. The vertices are enumerated starting from those located at the interior of S0 (with corresponding position vectors v1,...,vNI) followed by those located at S0 (with corresponding position vectors vNI+1,...,vNI+NB). The vector V3(NI+NB) is formed by concatenating the vertex position vectors vj,j=1,...,NI+NB, as follows: 
V=[v1vNI+NB]
(2)
The line segment coincident with a creased fold in S0, simply denoted as fold line, is defined by its start-point and end-point (both points corresponding to vertices by Definition 2.3). As such, to define the fold lines in S0, the matrix Ĉ{1,0,1}NF×(NI+NB) with elements Ĉij is introduced 
Ĉij={1;vjisthepositionoftheithfoldlinestart-point1;vjisthepositionoftheithfoldlineend-point0;otherwise,i.e.,vjisnotconnectedtotheithfoldline
(3)

Let nj, j=1,...,NI, be the number of fold lines incident to each interior vertex (i.e., those vertices located at the interior of S0). Also, let mjkspan(e1,e2),j=1,...,NI,k=1,...,nj, be the vector along the length of the kth fold line incident to the jth interior vertex that emanates from such a vertex (see Fig. 4).

Definition 2.4.The angle between the vectoryspan(e1,e2)ande1,which starts ate1and is measured in the counterclockwise direction, is denotedφ(y):span(e1,e2)[0,2π)and is determined as follows:2 
φ(y):={cos1(yy·e1);e2·y02πcos1(yy·e1);e2·y<0
(4)
For each interior vertex, its associated vectors mjk,k=1,...,nj, are arranged in a counterclockwise ordering such that φ(mj1)<φ(mj1)<<φ(mjnj)   j{1,...,NI}. Let Mj3nj be the vector constructed by concatenating the vectors mjk,k=1,...,nj, as follows: 
Mj=[mj1mjnj]
(5)
The matrices Cj{1,0,1}nj×NF,j=1,...,NI, with elements Ckij are used for the identification and ordering of the folds incident to the jth interior vertex and are defined as follows: 
Ckij={1;mjkisavectoralongtheithfoldlineandemanatesfromtheithfoldlinestart-point1;mjkisavectoralongtheithfoldlineandemanatesfromtheithfoldlineend-point0;otherwise 
(6)
The mapping from the vertex position vectors vj,j=1,...,NI+NB, to the vectors mjk,k=1,...,nj, is then compactly given as follows:3 
Mj=((CjĈ)I3)V
(7)

where In denotes the n×n identity matrix.

The face corner angles surrounding each interior vertex are denoted as αjk(0,2π],j=1,...,NI,k=1,...,nj, and are calculated as follows (see Fig. 4): 
αjk={φ(mjk+1)φ(mjk);k=1,...,nj12π+φ(mj1)φ(mjk);k=nj
(8)

Constraints.

Once the geometry of the fold pattern is defined, constraints on the fold shape variables (corresponding to the fold angles for creased folds) are formulated such that every current configuration is valid according to Definition 2.1.4 In addition to constraints for valid configurations, developability [40] is also conventionally imposed in origami. Developability allows a surface to be flattened onto a plane without stretching or overlapping. A developable surface has zero Gaussian curvature everywhere [79]. Since valid configurations of the sheet consist of planar faces joined at straight creased folds, the only location where the Gaussian curvature is nontrivially zero is at the interior fold intersections. The definition of Gaussian curvature as the product of the two principal curvatures is not valid at a singular interior fold intersection [80]; hence, the discrete Gaussian curvature, denoted Kj, is considered [54,80,81]. It is given as 2π less the sum of face corner angles surrounding each interior fold intersection [82], and it must be zero for the fold intersection to be developable 
Kj=2πk=1njαjk=0      j{1,...,NI}
(9)

In this work, the angles αjk are defined in S0, which is planar and having its faces only overlapping at their shared boundary edges, and thus, these angles clearly sum to 2π for each interior vertex. Therefore, the developability constraint in Eq. (9) is satisfied in S0. No further consideration of this constraint is required since the face corner angles are constant during the deformation history of the sheet (since the faces undergo only rigid deformations for valid configurations), and thus, they hold their associated values αjk as defined in S0. As a consequence, developability is assured for any valid current configuration St.

Nontrivial and important are the constraints on the fold angles that define the constrained configuration space. The set of constraints for θ̂i,i=1,...,NF, required for a valid configuration can be formulated as a set of constraints for the fold angles associated with the folds incident to each interior vertex (assuming that the sheet has no holes) [32].

The fold angle associated with the kth fold incident to the jth interior vertex is denoted θjk. The vectors θjnj,j=1,...,NI, are formed by collecting the fold angles θjk,k=1,...,nj, as follows: 
θj=[θj1θjnj]
(10)
and the mapping from the vector θ̂ with elements corresponding to the fold angles of the sheet (defined in Eq. (1)) to each vector θj is given as follows: 
θj=|Cj|*θ̂
(11)

where the elements of the matrix Cj are defined in Eq. (6), and |·|*:m×n0m×n denotes the elementwise absolute value of a matrix where [|Y|*]ij=|Yij|.

Let R1(ϕ)3×3 be the transformation matrix associated with a rotation by ϕ about an axis of rotation aligned to e1 
R1(ϕ):=[1000cos(ϕ)sin(ϕ)0sin(ϕ)cos(ϕ)]
(12)
and R3(ϕ)3×3 be the transformation matrix associated with a rotation by ϕ about an axis of rotation aligned to e3 
R3(ϕ):=[cos(ϕ)sin(ϕ)0sin(ϕ)cos(ϕ)0001]
(13)

The fold angles associated with the folds incident to each interior vertex θjk,k=1,...,nj, must be constrained to prevent tearing of closed strip of faces joined to each interior vertex [36,37]. For this purpose, the following constraint is proposed:

Proposition 2.1.For the initially closed strip of faces joined to the jth interior vertex to remain closed with each face undergoing a rigid deformation, the following constraint must hold [ib33,ib3232 ]: 
Rj:=k=1njR1(θjk)R3(αjk)=I3
(14)

It is shown in Appendix  A that this constraint for origami with creased folds is a special case of the constraints for origami with smooth folds to be presented in Sec. 3.3.

Corollary 2.1.

  • (i)

    If the jth interior vertex has a single incident creased fold, it allows for a valid configuration ifθj1=0.

  • (ii)

    If the jth interior vertex has two incident creased folds, it allows for a valid configuration if (ii,1)αj1π,θj1=θj2=0or (ii,2)αj1=π,θj1=θj2.

  • (iii)

    If the jth interior vertex has three incident creased folds, it allows for a valid configuration if (iii,1)αj1π,αj2π,αj3π,θj1=θj2=θj3=0, or (iii,2)αj1=π,θj1=θj2,θj3=0, or (iii,3)αj2=π,θj2=θj3,θj1=0, or (iii,4)αj3=π,θj3=θj1,θj2=0.

Proof. Case (i): The matrix Rj associated with an interior vertex having one incident creased fold is Rj=R1(θj1)R3(αj1). The angle αj1 is equal to 2π for a single fold incident to an interior vertex, thus R3(αj1)=I3. This requires R1(θj1)=I3 which in the domain θj1[π,π] holds true only when θj1=0. Cases (ii) and (iii) can be verified in a similar manner.□

Thus, for nontrivial folding deformation, any interior vertex must have at least four incident creased folds. It can be shown that each face in the sheet undergoes a rigid deformation and no tearing occurs provided the constraint in Eq. (14) is satisfied for each interior vertex of the sheet (see Refs. [36,37,83] for details). Therefore, the constraints on the fold angles defining the constrained configuration space for a sheet with creased folds (excluding self-intersection avoidance) are then the following: 
Rj=I3    j{1,...,NI}
(15)

For the sake of brevity, the formulation of the map from the reference configuration to a valid current configuration, constructed using the set of the shape variables of all the folds in the sheet, is not included here for origami creased folds and can be found in Refs. [36,37]. It is reminded that the model for origami with creased folds reviewed in this section is a special case of the more general model for origami with smooth folds presented in Sec. 3. Thus, the map from reference to current configurations for origami with smooth folds provided in Sec. 3.4 covers such a map for origami with creased folds.

Origami With Smooth Folds

The model for origami presented in the preceding section is based on the assumption of creased folds. Such a simplification may not be appropriate for structures having folds of non-negligible thickness or constructed from materials that do not provide sufficient strain magnitudes to generate the high curvatures needed for a creased idealization. For such structures, the folded regions may not be accurately represented as creases but instead as bent sheet regions exhibiting higher-order geometric continuity (termed smooth folds in this work). A novel model for the kinematic response of origami structures having smooth folds is presented in this section.

The approach used to develop a model for origami with smooth folds follows that outlined at the beginning of Sec. 2. First, the sheet and the shape variables associated with the smooth folds in the sheet are described. The layout of the smooth folds in the sheet (i.e., the fold pattern) is determined by vertices (start-points and end-points of the smooth fold centerlines in a planar reference configuration), their connectivity, and the initial width of each smooth fold. Subsequently, constraints on the fold shape variables that define the constrained configuration space are derived. The continuous motion of the sheet represents a continuous path in such a constrained configuration space.

The studied continuum body is denoted as the sheet which is a three-dimensional, orientable, path-connected surface with boundary (same properties of a sheet in origami with creased folds, see Sec. 2). For origami with smooth folds, the sheet is divided into various surface subdomains denoted as faces, smooth folds, and fold intersections. Every surface subdomain comprising the sheet has the same aforementioned characteristics of the sheet (i.e., they are three-dimensional, orientable, path-connected surfaces with boundary). Following the notation of Sec. 2, the orthonormal vectors ei3,i=1,2,3, with e3:=e1×e2 form the basis {e1,e2,e3} that defines the fixed global coordinate system.

The reference configuration of the sheet is denoted S0 and is defined such that it is contained in the plane spanned by e1 and e2 with its surface subdomains not overlapping each other, except at their shared boundary edges. The configuration of the NP faces, NF smooth folds, and NI fold intersections subdomains in S0 is denoted P0i,i=1,...,NP,F0i,i=1,...,NF, and I0i,i=1,...,NI, respectively. Therefore, S0=(i=1NPP0i)(i=1NFF0i)(i=1NII0i). The side of S0 with normal e3 is selected as the positive side of the sheet. An example of a sheet with smooth folds in its reference configuration S0 is shown in Fig. 5(a).

A current configuration of the sheet is denoted St where t is a parameter that indicates the history of deformation from the reference configuration (t = 0) to a current configuration (t > 0). Refer to Fig. 5(b) for an example. The configuration of the NP faces, NF smooth folds, and NI fold intersections comprising St is denoted Pti,i=1,...,NP,Fti,i=1,...,NF, and Iti,i=1,...,NI, respectively. Thus, St=(i=1NPPti)(i=1NFFti)(i=1NIIti). As stated in Sec. 1, the focus of this work is on rigid origami-based structures. As such, the following definition provides the characteristics of a valid current configuration for origami with smooth folds:

Definition 3.1.Valid configuration:A valid current configurationSthas the following characteristics: (i) the faces undergo only rigid deformations, (ii) the sheet is not torn (initially joined surface subdomains of the sheet remain joined), and (iii) the sheet does not self-intersect.

The characteristics of a valid configuration provided in Definition 3.1 are the same as those of a valid configuration for a sheet with creased folds presented in Definition 2.1. However, it is noted that unlike a sheet with creased folds, a sheet with smooth folds is also comprised of other surface subdomains besides the faces (smooth folds and fold intersections) for which bending and stretching are permitted.

Definition 3.2.Smooth folds:The smooth folds are ruled surfaces5 of the following form: 
Fti(ζ1,ζ2)=cti(ζ1)+ζ2hti,dcti(ζ1)dζ1·hti=0
(16)

whereFti(ζ1,ζ2)3is a parameterization ofFti. Without loss of generality, the parameters ζ1 and ζ2 are contained in the intervals[1,1]and[0,1], respectively.

An example of a smooth fold is shown in Fig. 6. In Definition 3.2, hti3 provides the direction of the rulings comprising Fti while cti(ζ1)3 is the parametric ruled surface directrix curve that defines the cross section of Fti. The curve parameterized by cti(ζ1) is contained in a plane orthogonal to hti as stated in Eq. (16). It is assumed that hti is constant for all configurations. As a consequence, the only nonrigid deformations allowed for the smooth folds are achieved through continuous bending or stretching of its cross section defined by cti(ζ1). To simplify the notation, the dependence of Fti(ζ1,ζ2),cti(ζ1), and hti on t is left implicit for the remainder of the paper and these vectors are denoted as Fi(ζ1,ζ2),ci(ζ1), and hi, respectively.

The reference configuration of the smooth folds is simplified here to a rectangular shape and their deformation only includes stretching and bending of the fold cross section. However, the modeling framework proposed in this work can be extended for the consideration of folds having a trapezoidal reference configuration (that deform into conical sections) or folds that exhibit torsional deformation. Such extensions are recommended for future studies.

Remark 3.1. Deformation of a sheet with creased folds is the special limiting case of the more general deformation of a sheet with smooth folds. Specifically, creased folds are obtained when the curve parameterized by ci(ζ1) is degenerated to a single point, thereby degenerating the smooth fold surface Fti to a single straight line segment.

Each smooth fold is joined to a face at each of its two boundary rulings (i.e., Fi(1,ζ2) and Fi(1,ζ2), refer to Fig. 6). The remaining boundaries of Fti (i.e., Fi(ζ1,0) and Fi(ζ1,1)) are either at the boundary of St or joined to a fold intersection. A parameterization for the fold intersections (Iti,i=1,...,NI) is not provided here. It is noted, however, that the model proposed herein is independent of the parameterization of Iti.

A nonrigid body deformation of the sheet is achieved by rotating pairs of faces joined to smooth folds relative to one another in such a manner that the sheet only attains valid configurations during such a deformation. One of the shape variables associated with a smooth fold describes the relative rotation between the two faces joined by such a fold and is denoted as fold angle.

Definition 3.3.Fold angle:The fold angleθ̂i(t)is defined as π minus the dihedral angle between the positive sides of the two faces joined to the ith smooth fold.

The dependence of the fold angles on t is left implicit for the remainder of the paper. The vector θ̂NF was constructed by collecting the fold angles θ̂i,i=1,...,NF, and is defined in Eq. (1). Schematics showing the concept of fold angle for smooth folds are provided in Fig. 6.

Geometry of Smooth Folds.

This section presents the detailed geometrical description of smooth folds. The conditions required for various orders of geometric continuity and particular formulations for ci(ζ1) are also provided.

The distance between the two end-points of ci(ζ1) with position vectors ci(1) and ci(1) in any configuration is denoted ŵi 
ŵi:=ci(1)ci(1)
(17)

A schematic of a smooth fold cross section showing ŵi is provided in Fig. 7. The vector ŵNF is constructed by collecting the variables ŵi as follows:

 
ŵ=[ŵ1ŵNF]
(18)
The fold widthŵi0 is the value of ŵi at t = 0. The vector ŵ0NF is constructed as follows: 
ŵ0=[ŵ10ŵNF0]
(19)

The fold-attached orthonormal vectors êji3,i=1,...,NF, j = 1, 2, 3, with ê3i:=ê1i×ê2i form the bases {ê1i,ê2i,ê3i} that define the local fold coordinate system of each smooth fold Fti. The origin of this coordinate system is located at (1/2)(ci(1)+ci(1)). The director vector hi is aligned to ê1i (i.e., hi·ê1i=hi) and (ci(1)ci(1))span(ê2i).

The face adjacent to a smooth fold at the boundary Fi(1,ζ2) makes an angle of âiθ̂i with ê2i in the plane spanned by ê2i and ê3i, and the face adjacent to a smooth fold at the boundary Fi(1,ζ2) makes an angle of (1âi)θ̂i with ê2i in such a plane. This is shown schematically in Fig. 7. The vector âNF is constructed by collecting the variables âi as follows: 
â=[â1âNF]
(20)
Let ĉi(ζ1) be the parametric curve ci(ζ1) expressed in the fold coordinate system of Fti 
ĉi(ζ1)=[ê1iê2iê3i][ci(ζ1)12(ci(1)+ci(1))]
(21)
If ĉi(ζ1) is at least first-order differentiable for ζ1[1,1], the total arc length of the fold cross section, denoted sitot, is determined as follows: 
sitot=11dĉi(ζ1)dζ1dζ1
(22)
and the arc length si(ζ1) measured from ζ1=1 to any value of ζ1[1,1] is determined as follows: 
si(ζ1)=1ζ1dĉi(ζ)dζdζ
(23)
The tangent direction of the parametric curve ĉi(ζ1) is determined by the unit tangent vector ti(ζ1)span(ê2i,ê3i) and is defined at the points where ĉi(ζ1) is at least first-order differentiable 
ti(ζ1)=dĉi(ζ1)dζ1/dĉi(ζ1)dζ1
(24)
The curvature κi(ζ1) is defined at the points where ĉi(ζ1) is at least second-order differentiable 
κi(ζ1)=dĉi(ζ1)dζ1×d2ĉi(ζ1)dζ12/dĉi(ζ1)dζ13
(25)
and the signed curvature is given as follows:6 
κis(ζ1)=κi(ζ1)sgn((dĉi(ζ1)dζ1×d2ĉi(ζ1)dζ12)·ê1i)
(26)
The order of geometric continuity of St at the interior rulings of Fti is determined by the order of geometric continuity of ĉi(ζ1),ζ1(1,1), while that at the joints with the planar faces adjacent to Fti depends on the particular values of ĉi(ζ1) and its derivatives at ζ1=±1. For instance, G0 continuous joints with the faces adjacent to Fti require the following conditions on ĉi(ζ1) at ζ1=±1 (refer to Fig. 7): 
ĉi(1)=[0ŵi20]=:ĉL0i,ĉi(1)=[0ŵi20]=:ĉR0i
(27)
Continuity of the unit tangent vector ti(ζ1) at ζ1=±1 is required for G1 continuous joints with the planar faces adjacent to Fti in addition to G0 continuity [84,85]. The following values of ti(ζ1) at ζ1=±1 are then required for G1 continuity in addition to those conditions of Eq. (27): 
ti(1)=[0cos(âiθ̂i)sin(âiθ̂i)],  ti(1)=[0cos((1âi)θ̂i)sin((1âi)θ̂i)]
(28)
Therefore, the following conditions on the first derivatives of ĉi(ζ1) at ζ1=±1 are required for G1 continuity at the joints with the planar faces adjacent to Fti : 
dĉi(ζ1)dζ1|ζ1=1=βL1i[0cos(âiθ̂i)sin(âiθ̂i)]=:ĉL1idĉi(ζ1)dζ1|ζ1=1=βR1i[0cos((1âi)θ̂i)sin((1âi)θ̂i)]=:ĉR1i
(29)

where βL1i,βR1i>0.

Continuity of the curvature vector (or the signed curvature) is required for G2 continuity [85] in addition to G1 continuity. This requires the curvature of ĉi(ζ1) at ζ1=±1 to be zero since Fti is connected to planar faces at its boundary rulings. Zero curvature at the boundary rulings of Fti requires the following according to Eq. (25): 
(dĉi(ζ1)dζ1×d2ĉi(ζ1)dζ12)|ζ1=1=03(dĉi(ζ1)dζ1×d2ĉi(ζ1)dζ12)|ζ1=1=03
(30)
where 0n is the zero vector in n. Considering Eqs. (29) and (30), the following conditions on the second derivatives of ĉi(ζ1) are needed in addition to the conditions provided in Eqs. (27) and (29) for G2 continuity at the joints with the planar faces adjacent to Fti : 
d2ĉi(ζ1)dζ12|ζ1=1=βL2i[0cos(âiθ̂i)sin(âiθ̂i)]=:ĉL2id2ĉi(ζ1)dζ12|ζ1=1=βR2i[0cos((1âi)θ̂i)sin((1âi)θ̂i)]=:ĉR2i
(31)

where βL2i,βR2i. Conditions on higher-order derivatives of ĉi(ζ1) required for higher-order geometric continuity can be provided in a similar manner [84,85].

Fold Shape Examples.

Polynomials of the minimum order required to satisfy the previous conditions for continuity of ĉi(ζ1) and its derivatives at ζ1=±1 are used to define this parametric curve. Hermite interpolation polynomials [86] are used to represent ĉi(ζ1) in this work. Alternative representations (e.g., Bezier curves) are also applicable as long as they satisfy the conditions required for the considered order of geometric continuity. For G1 continuity at the boundary rulings of Fti,ĉi(ζ1) is expressed as follows: 
ĉi(ζ1)=h30(ζ1)ĉL0i+h31(ζ1)ĉR0i+h32(ζ1)ĉL1i+h33(ζ1)ĉR1i
(32)
where the utilized cubic Hermite interpolation polynomials h3i(ζ),i=0,...,3, are given as follows: 
h30(ζ)=14ζ334ζ+12h31(ζ)=14ζ3+34ζ+12h32(ζ)=14ζ314ζ214ζ+14h33(ζ)=14ζ3+14ζ214ζ14
(33)
For G2 continuity at the boundary rulings of Fti,ĉi(ζ1) is expressed as follows: 
ĉi(ζ1)=h50(ζ1)ĉL0i+h51(ζ1)ĉR0i+h52(ζ1)ĉL1i +h53(ζ1)ĉR1i+h54(ζ1)ĉL2i+h55(ζ1)ĉR2i
(34)
where the utilized quintic Hermite interpolation polynomials h5i(ζ),i=0,...,5, are given as follows: 
h50(ζ)=316ζ5+58ζ31516ζ+12h51(ζ)=316ζ558ζ3+1516ζ+12h52(ζ)=316ζ5+116ζ4+58ζ338ζ2716ζ+516h53(ζ)=316ζ5116ζ4+58ζ3+38ζ2716ζ516h54(ζ)=116ζ5+116ζ4+18ζ318ζ2116ζ+116h55(ζ)=116ζ5+116ζ418ζ318ζ2+116ζ+116
(35)

Fold Pattern.

The layout of the smooth folds in S0 is presented in a fold pattern. To describe the fold pattern, the vertices are first introduced.

Definition 3.4.Vertices:The vertices are the start-points and end-points of the line segments coincident with the centerlines of the smooth folds inS0. Each vertex has an associated position vector denotedvjspan(e1,e2).

As in Sec. 2.1, the number of vertices located at the interior of S0 is denoted NI, and the number of vertices located at S0 or outside S0 is denoted NB. The vertices are enumerated starting from those located at the interior of S0 (with corresponding position vectors v1,...,vNI) followed by those at the located at S0 or outside S0 (with corresponding position vectors vNI+1,...,vNI+NB). The vertices in origami with smooth folds are used only to indicate end-points of the fold centerlines in S0. The points of the sheet coincident with the vertices in S0 do not occupy any especial location in St, t > 0. This is in contrast to origami with creased folds where the points of the sheet coincident with the vertices in S0, which correspond to fold line end-points, and also correspond to fold line end-points in St, t > 0 (see Sec. 2).

Following Sec. 2.1, the vector V3(NI+NB) formed by concatenating the vertex position vectors is provided in Eq. (2). The definition of the matrix Ĉ{1,0,1}NF×(NI+NB) that identifies which vertices are the start-points and end-points of the fold centerlines follows that provided in Eq. (3). Let v̂i1,v̂i2span(e1,e2),i=1,...,NF, be the position vectors of the vertices from which each fold centerline emanates and ends, respectively. The vectors V̂1,V̂23NF are constructed by concatenating the vectors v̂i1,v̂i2,i=1,...,NF, as follows: 
V̂1=[v̂11v̂NF1],    V̂2=[v̂12v̂NF2]
(36)
The mappings from the collection of all the vertex position vectors V to those corresponding to the position vectors of the start-points and end-points of the fold centerlines (V̂1 and V̂2, respectively) are given as follows: 
V̂1=(12(|Ĉ|*Ĉ)I3)VV̂2=(12(|Ĉ|*+Ĉ)I3)V
(37)
The four corner points of F0i having associated position vectors p̂1i,p̂2i,p̂3i,p̂4ispan(e1,e2),i=1,...,NF, are then determined as follows: 
p̂1i=v̂i1ŵi02(e3×v̂i2v̂i1v̂i2v̂i1)+r̂1iv̂i2v̂i1v̂i2v̂i1p̂2i=v̂i2ŵi02(e3×v̂i2v̂i1v̂i2v̂i1)r̂2iv̂i2v̂i1v̂i2v̂i1p̂3i=v̂i2+ŵi02(e3×v̂i2v̂i1v̂i2v̂i1)r̂2iv̂i2v̂i1v̂i2v̂i1p̂4i=v̂i1+ŵi02(e3×v̂i2v̂i1v̂i2v̂i1)+r̂1iv̂i2v̂i1v̂i2v̂i1
(38)
where r̂1i,r̂2i and the resulting fold length orthogonal to the width direction must be positive 
v̂i2v̂i1r̂1ir̂2i>0       i{1,...,NF}
(39)

The geometric parameters that define the corner points of F0i are shown in Fig. 8.

The vectors P̂j3NF,j=1,...,4, are constructed by concatenating the vectors p̂ji,i=1,...,NF, as follows: 
P̂j=[p̂j1p̂jNF]
(40)

As in Sec. 2.1, let nj, j=1,...,NI, be the number of fold centerlines incident to each interior vertex (i.e., those vertices located at the interior of S0). Also, let mjkspan(e1,e2),j=1,...,NI,k=1,...,nj, be the vector along the length of the kth fold centerline incident to the jth interior vertex that emanates from this vertex. The vectors mjk,k=1,...,nj, have a counterclockwise ordering, i.e., φ(mj1)<φ(mj1)<<φ(mjnj)j{1,...,NI}, where φ(·) is defined in Eq. (4).

The mapping from the vertex position vectors vj,j=1,...,NI+NB, to the vectors mjk,k=1,...,nj, is provided in Eq. (7). The angles between adjacent fold centerlines intersecting at a common interior vertex αjk,j=1,...,NI,k=1,...,nj, are calculated using Eq. (8). A schematic showing the vectors mjk and the angles αjk is provided in Fig. 9.

Constraints.

As in conventional origami with creased folds [36,37], constraints are required for origami with smooth folds to ensure that every current configuration St is valid7 (according to Definition 3.1).

The angles αjk, k=1,...,nj, satisfy the following constraint since they are defined in S0 : 
Kj:=2πk=1njαjk=0     j{1,...,NI}
(41)

Since the faces undergo only rigid deformations, the angles αjk are constant during the deformation history of the sheet, i.e., independent of t, and therefore, Eq. (41) holds at every configuration St. Equation (41) represents the developability constraint for origami with creased folds (see Eq. (9)). However, pointwise isometric deformation is relaxed here for origami with smooth folds since stretching is permitted within the subdomains Fti and Iti, and thus, only S0 is in general pointwise developable. It is remarked that isometry is assumed for the smooth folds Fti in the direction of hi (see Fig. 6). However, the arc length of the curve parameterized by ĉi(ζ1) may change in general during deformation, and thus, stretching of the smooth folds in such a direction may be allowed. Refer to Definition 3.2 and its subsequent discussion.

The variables describing the deformation associated with the folding of the smooth folds are constrained such that every configuration St is valid. These variables correspond to the fold angle θ̂i, the distance ŵi between the end-points of the cross section parametric curve ĉi(ζ1), and âi,i=1,...,NF (refer to Sec. 3.1). Following Sec. 2.2, a new set of constraints for the variables associated with the smooth folds adjacent to each interior fold intersection are now proposed.

The variables θjk, wjk, and ajk are those associated with the kth smooth fold adjacent to I0j. The vectors θjnj,j=1,...,NI, constructed by collecting the fold angles θjk, k=1,...,nj, are defined in Eq. (10). The vectors wj,ajnj,j=1,...,NI, are constructed by collecting the variables wjk and ajk, k=1,...,nj, and are defined as follows: 
wj=[wj1wjnj],aj=[aj1ajnj]
(42)
The mapping from the vector θ̂ with elements corresponding to the fold angles of the sheet to each vector θj with elements corresponding to the fold angles of the smooth folds adjacent to I0j is provided in Eq. (11). The mapping from the vector ŵ, defined in Eq. (18), to each vector wj is the following: 
wj=|Cj|*ŵ
(43)
where the definition of the matrix Cj in Eq. (6) holds in the context of connectivity of fold centerlines. Taking into account the orientation of the adjacent smooth folds with respect to the considered fold intersection I0j (i.e., whether the interior vertex associated with I0j is the start-point or the end-point of the adjacent smooth fold centerline), the mapping from the vector â, defined in Eq. (20), to each vector aj is as follows: 
aj=[CjACj][â1]
(44)
where ACjnj is a vector with elements ACkj that are determined as follows: 
ACkj=12(1i=1NFCkij)
(45)

Note that if the vector mjk is coincident with and has the same orientation as the ith fold centerline, then ajk=âi. Conversely, if the vector mjk is coincident with and has the opposite orientation as the ith fold centerline, then ajk=1âi. This is obtained from the mapping provided in Eq. (44).

The fold widths associated with the smooth folds adjacent to I0j are denoted wjk0,k=1,...,nj, and are required in the subsequent derivations. Let w0,jnj be the vector constructed by collecting the fold widths wjk0,k=1,...,nj. Such a vector and its mapping from ŵ0 (defined in Eq. (19)) are as follows: 
w0,j=[wj10wjnj0]=|Cj|*ŵ0
(46)

Let γj(η):[0,1]S0 be any simple closed path (i.e., γj(0)=γj(1)) enclosing I0j and crossing each smooth fold adjacent to I0j once. An example of a path γj(η) is shown in Fig. 10. The point having position γj(0)=γj(1) is defined such that it is located at the face adjacent to the smooth folds with corresponding vectors mj1 and mjnj. Also, the path γj(η) is defined such that it crosses the smooth folds with associated vectors mjk in counterclockwise order (i.e., mj1,mj2,...,mjnj).

The position vectors of the points where γj(η) crosses each boundary ruling of the smooth folds with associated vector mjk are denoted bLjkspan(e1,e2) (point where γj(η) enters the smooth fold) and bRjkspan(e1,e2) (point where γj(η) exits the smooth fold). This is shown schematically in Fig. 10(a).

Let Q1(ϕ)4×4 be the transformation matrix in homogeneous coordinates associated with a rotation by ϕ about an axis of rotation aligned to e1 
Q1(ϕ):=[1 00 00 cos(ϕ)sin(ϕ) 00 sin(ϕ)cos(ϕ) 00 00 1]=[R1(ϕ)03031]
(47)
and Q3(ϕ)4×4 be the transformation matrix associated with a rotation by ϕ about an axis of rotation aligned to e3 
Q3(ϕ):=[cos(ϕ)sin(ϕ) 0  0sin(ϕ)cos(ϕ) 0  000 1  000 0  1]=[R3(ϕ)03031]
(48)
Let T(b)4×4 be the matrix representing the transformation associated with a translation by vector b3 with elements bi 
T(b):=[1 0 0 b10 1 0 b20 0 1 b30 0 0 1]=[I3b031]
(49)
Considering an axis with direction aligned to a vector yspan(e1,e2) that crosses a point with position vector bspan(e1,e2), the transformation associated with a rotation by ϕ about such an axis can be represented as follows [36,37]: 
T(b)Q3(φ(y))Q1(ϕ)Q31(φ(y))T1(b)
(50)
Theorem 3.1.The transformation matrixLjkdescribing the deformation associated with the folding of the kth smooth fold crossed by the pathγj(η)is determined as follows:8 
Ljk=(T(bLjkgjk1)Q3(φ(mjk))Q1(ajkθjk)Q31(φ(mjk))T1(bLjkgjk1)T(bRjkgjk)Q3(φ(mjk))Q1((1ajk)θjk)Q31(φ(mjk))T1(bRjkgjk))
(51)
where the vectorsgjkspan(e1,e2),k=0,...,nj, account for the change in the distance between the boundary rulings of the smooth folds in a current configuration and are determined recursively as follows: 
gj0=03,gjk=gjk1+(wjk0wjk)(e3×mjkmjk)   k{1,...,nj}
(52)

Proof. Referring to Figs. 11(c) and 11(d), the transformation associated with the folding of the kth smooth fold crossed by γj(η) can be discretized by two consecutive rotation transformations. The first transformation corresponds to a rotation by (1ajk)θjk about an axis aligned to mjk and crossing a point with position vector bRjkgjk. Using Eq. (50), the transformation matrix associated with this rotation is the following:

 
(T(bRjkgjk)Q3(φ(mjk))Q1((1ajk)θjk)     Q31(φ(mjk))T1(bRjkgjk))
(53)
The second transformation corresponds to a rotation by ajkθjk about an axis aligned to mjk and a crossing point with position vector bLjkgjk1. The transformation matrix associated with this rotation has the form 
(T(bLjkgjk1)Q3(φ(mjk))Q1(ajkθjk)     Q31(φ(mjk))T1(bLjkgjk1))
(54)

The composition of the transformations shown in Eqs. (53) and (54) results in the transformation matrix Ljk provided in Eq. (51).

Referring to Figs. 11(a) and 11(b), the change in the distance between the boundary rulings of a smooth fold in a current configuration is given by (wjk0wjk). To define the position of the axes of rotation for the transformations composing Ljk, the vector l=1k1(wjl0wjl)[e3×(mjl/mjl)] must be subtracted from bLjk while l=1k(wjl0wjl)[e3×(mjl/mjl)] must be subtracted from bRjk. The recursive definition of the vectors gjk,k=0,...,nj, in Eq. (52) allows for a simplified form of these vector subtractions.□

Utilizing the block matrix expressions for the rotation and translation matrices given in Eqs. (47)(49), the transformation matrices Ljk provided in Eq. (51) can be partitioned into four blocks: ([Ljk]11block3×3, [Ljk]12block3×1, [Ljk]21block1×3, and [Ljk]22block1×1) as follows: 
[Ljk]11block=R3(φ(mjk))R1(θjk)R31(φ(mjk))
(55)
 
[Ljk]12block=bLjkgjk1+(R3(φ(mjk))R1(ajkθjk)R31(φ(mjk))(bRjkgjkbLjk+gjk1))(R3(φ(mjk))R1(θjk)R31(φ(mjk))(bRjkgjk))
(56)
 
[Ljk]21block=03
(57)
 
[Ljk]22block=1
(58)
The vectors wjk, ljkspan(e1,e2) are defined as follows: 
wjk:=bRjkgjkbLjk+gjk1
(59)
 
ljk:={bLjk+1bRjk;k=1,...,nj1bLj1bRjk;k=nj
(60)
The vectors wjk and ljk for the example in Fig. 10(a) are shown in Fig. 10(b). Let the vectors w̃jk and l̃jk be wjk and ljk, respectively, expressed in a coordinate system with the one-axis aligned to mjk and the three-axis aligned to e3 
w̃jk=R31(φ(mjk))wjk
(61)
 
l̃jk=R31(φ(mjk))ljk
(62)
Following the process presented in Refs. [36,37], to formulate constraints for the shape variables associated with the folds, the map from the reference to current configurations considering only the faces and smooth folds adjacent to an interior fold intersection is first constructed. The face containing the point with position γj(0) is assumed fixed in space (not translating or rotating) for the derivation of constraints. Let Xspan(e1,e2) be the position vector of a point in a face adjacent to I0j in the reference configuration and let x3 be the position vector of such a point in a current configuration. The map Xx is constructed as the composition of transformations Ljk associated with the folds crossed by the segment of path γj(η) that connects γj(0) to the face containing the point with initial position X 
[x1]=(k=1nyLjk)T1(gjny)[X1] =(k=1nyHjk)[X1]
(63)
where the matrices Hjk,j=1,...,NI,k=1,...,nj, are expressed as follows: 
Hjk=(T(bLjk)Q3(φ(mjk))Q1(ajkθjk)Q31(φ(mjk))T1(bLjk)T1((wjk0wjk)(e3×mjkmjk))T(bRjk)Q3(φ(mjk))Q1((1ajk)θjk)Q31(φ(mjk))T1(bRjk))
(64)

and ny is the number of smooth folds crossed by the segment of the path γj(η) that connects γj(0) and the face containing the point with position vector X. Note that x is the position vector of such a point in a current configuration determined by θjk, ajk, and wjk, k=1,...,nj. Since such a mapping is a composition of translation and rotation matrices, each face undergoes a rigid deformation (required for a valid configuration). In order to prevent tearing among the surface subdomains joined to Itj, the following constraints are proposed:

Theorem 3.2.For the initially closed strip of faces and smooth folds joined toI0jto remain closed with each face undergoing a rigid deformation, two constraints must hold. These are the rotation constraint 
Rj:=k=1njR1(θjk)R3(αjk)=I3
(65)
and the translation constraint 
dj:=k=1nj((l=1k1R1(θjl)R3(αjl))R1(ajkθjk)w̃jk+(l=1k1R1(θjl)R3(αjl))R1(θjk)l̃jk)=03
(66)
Proof. The deformation map of a point in the face containing the point with position γj(0) must be the identity transformation since such a face is assumed fixed (i.e., if ny = nj in Eq. (63) is considered, then x=X). This requires the following: 
I4=(k=1njHjk)=(k=1njLjk)T1(gjnj)
(67)
Utilizing Eqs. (55)(58), Eq. (67) can be partitioned into four blocks. The 11 block is the following: 
I3=j=1njR3(φ(mjk))R1(θij)R31(φ(mjk))=(R3(φ(mj1))(k=1nj1R1(θjk)R3(αjk))R1(θjnj)R31(φ(mjnj)))
(68)
where the following equality was used: 
R3(αjk)={R31(φ(mjk))R3(φ(mjk+1));k=1,...,nj1R31(φ(mjk))R3(φ(mj1));k=nj
(69)
Multiplying the last expression in Eq. (68) by R31(φ(mj1)) from the left and by R3(φ(mj1)) from the right, the following is obtained: 
I3=((k=1nj1R1(θjk)R3(αjk))R1(θjnj)R31(φ(mjnj))R3(φ(mj1)))=(k=1nj1R1(θjk)R3(αjk))R1(θjnj)R3(αjnj)=k=1njR1(θjk)R3(αjk)=Rj
(70)
cf. Eq. (65). The 12 block of Eq. (67) is the following: 
03=k=1nj((l=1k1R3(φ(mjl))R1(θjl)R31(φ(mjl)))(bLjkgjk1+(R3(φ(mjk))R1(ajkθjk)R31(φ(mjk))(bRjkgjkbLjk+gjk1))(R3(φ(mjk))R1(θjk)R31(φ(mjk))(bRjkgjk))))(k=1njR3(φ(mjk))R1(θjk)R31(φ(mjk)))gjnj
(71)
Substituting the second expression of Eq. (68) and Eqs. (61) and (62) into Eq. (71), the following is obtained: 
03=k=1nj((l=1k1R(φ(mjl))R1(θjl)R31(φ(mjl)))(R3(φ(mjk))R1(ajkθjk)w̃jk+R3(φ(mjk))R1(θjk)l̃jk))
(72)
Finally multiplying both sides of the previous expression by R31(φ(mj1)) and simplifying using Eq. (69), the following is obtained: 
03=k=1nj((l=1k1R1(θjl)R3(αjl))R1(ajkθjk)w̃jk+(l=1k1R1(θjl)R3(αjl))R1(θjk)l̃jk)=dj
(73)

cf. Eq. (66). The 21 and the 22 blocks of the right side of Eq. (67) are equal to 03 and 1, respectively.□

Corollary 3.1.

  • (i)

    IfI0jhas a single adjacent smooth fold, it allows for a valid configuration ifθj1=0.

  • (ii)

    IfI0jhas two adjacent smooth folds, it allows for a valid configuration if (ii,1)αj1π,θj1=θj2=0or (ii,2)αj1=π,θj1=θj2.

  • (iii)

    IfI0jhas three adjacent smooth folds, it allows for a valid configuration if (iii,1)αj1π,αj2π,αj3π,θj1=θj2=θj3=0, or (iii,2)αj1=π,θj1=θj2,θj3=0, or (iii,3)αj2=π,θj2=θj3,θj1=0, or (iii,4)αj3=π,θj3=θj1,θj2=0.

Since the rotation constraint for smooth folds (Rj=I3, cf. Eq. (65)) also holds for creased folds (cf. Eq. (14)), the proof of this corollary is the same as that provided for Corollary 2.1.

The selection of the vectors bLjk,bRjk,k=1,...,nj is not unique and depends on the particular choice of the closed path γj(η) (see Fig. 10(a)). The corner points of the smooth folds defined in Eq. (38) provide a simple choice for the points where γj(η) crosses the boundary rulings of each smooth fold adjacent to I0j. Thus, they can be used to define bLjk,bRjk,k=1,...,nj. First, let BLj,BRj3nj be the vectors constructed by concatenating the vectors bLjk,bRjk,k=1,...,nj, as follows: 
BLj=[bLj1bLjnj],BRj=[bRj1bRjnj]
(74)
Taking into account the orientation of the adjacent smooth folds with respect to the considered fold intersection I0j (i.e., whether the interior vertex associated with I0j is the start-point or the end-point of the adjacent smooth fold centerline), the mapping from the corner points of the smooth folds to BLj and BRj is compactly given as follows: 
[BLjBRj]=[12(|Cj|*+Cj)I312(|Cj|*Cj)I312(|Cj|*Cj)I312(|Cj|*+Cj)I3][P̂1P̂4]
(75)

where the vectors P̂i are defined in Eq. (40). Note that the mapping provided in Eq. (75) results in bLjk=p̂1i,bRjk=p̂4i if mjk is coincident with and has the same orientation as the ith smooth fold centerline and bLjk=p̂4i,bRjk=p̂1i if mjk is coincident with and has the opposite orientation as the ith smooth fold centerline.

Finally, the constraints on the fold shape variables that allow for a valid configuration (excluding self-intersection avoidance) for a sheet with smooth folds are then the following: 
Rj=I3,dj=03j{1,...,NI}
(76)

cf. Eq. (15).

In general, the variables θ̂i,ŵi, and âi for each individual smooth fold represent its shape variables (i.e., degrees-of-freedom) that describe its configuration. If certain assumptions hold regarding the deformation of each individual smooth fold, the associated variables θ̂i,ŵi, and âi may not vary independently but rather relations among them can be derived. For simplicity in the implementation of the proposed model, assumptions on the extensibility and curvature field of the curve ĉi(ζ1) are taken such that the overall deformation of a smooth fold becomes a function of the fold angle θ̂i and the fold width ŵi0 (i.e., ŵi=ŵi(θ̂i,ŵi0),âi=âi(θ̂i,ŵi0) for each fold). Such a process involves nondimensionalization of the parametric curve ĉi(ζ1) and is presented in detail in Appendix  B for the case of smooth folds exhibiting G2 continuity (Eq. (34)).

Folding Map.

The formulation of the deformation map that relates the reference and current configurations (termed folding map here) is provided in this section. A folding map considering only the smooth folds and faces adjacent to an interior fold intersection was provided in Eq. (63). Here such a formulation is extended for the folding map of all the smooth folds and faces in the sheet.

First, an arbitrary face in the sheet is assumed fixed in its reference configuration.9 Let γj(η):[0,1]S0,j=1,...,NP, be paths connecting the fixed face to P0j,j=1,...,NP, respectively (see Fig. 12 for an example). The paths γj(η),j=1,...,NP, may not cross any fold intersection (i.e., they pass only through faces and smooth folds of S0). Following the formulation of Eq. (63), the map for points located in the faces is constructed as the composition of transformations associated with the smooth folds crossed by the path γj(η) between γj(0) (located at the fixed face) and γj(1) (located at P0j). Each path γj(η) crosses a number of pj smooth folds.

The path γj(η) crosses a smooth fold positively if it enters such a fold from the ruling Fi(1,ζ2) and exits at the ruling Fi(1,ζ2). Conversely, the path γj(η) crosses a smooth fold negatively if it crosses such a fold in the opposite direction. The matrices Cj{1,0,1}pj×NF,j=1,...,NP, with elements Ckij are used for the identification and ordering of the folds crossed by γj(η) 
Ckij={1;F0iisthekthfoldcrossedbyγj(η)andispositivelycrossed1;F0iisthekthfoldcrossedbyγj(η)andisnegativelycrossed0;otherwise
(77)
Let mjk be the vector along the length of the centerline of the kth smooth fold crossed by γj(η). Such a vector has the same orientation as the fold centerline if γj(η) crosses the fold positively and opposite orientation if γj(η) crosses the fold negatively. The vector Mj3pj is constructed by concatenating the vectors mjk,k=1,...,pj, as follows: 
Mj=[mj1mjpj]
(78)
The mapping from the vertex position vectors collected in the vector V to each vector Mj is compactly given as follows: 
Mj=((CjĈ)I3)V
(79)
Let θjk,wjk, and ajk,k=1,...,pj, be the shape variables associated with the ordered smooth folds crossed by γj(η). The vectors θj,wj,ajpj are constructed by collecting the variables θjk,wjk, and ajk,k=1,...,pj, as follows: 
θj=[θj1θjpj],wj=[wj1wjpj]aj=[aj1ajpj]
(80)
Following the formulation presented in Sec. 3.3, the vectors θj,wj, and aj are determined from the shape variables of the smooth folds in the sheet as follows (cf. Eqs. (11), (43), and (44)): 
θj=|Cj|*θ̂,wj=|Cj|*ŵ,aj=[CjACj][â1]
(81)
where ACjpj is a vector with elements ACkj determined as follows (cf. Eq. (45)): 
ACkj=1212i=1NFCkij
(82)
The fold widths associated with the smooth folds crossed by γj(η) are denoted wjk0,k=1,...,nj, and are required in the formulation of the folding map. As such, let w0,jpj be the vector constructed by collecting the fold widths wjk0,k=1,...,pj. Such a vector and its mapping from ŵ0 (defined in Eq. (19)) are as follows: 
w0,j=[wj10wjnj0]=|Cj|*ŵ0
(83)
The position vectors of the points where γj(η) crosses the boundary rulings of the smooth folds are denoted bLjk3 (point where γj(η) enters the smooth fold) and bRjk3 (point where γj(η) exits the smooth fold). Let BLj,BRj3pj be the vectors constructed by concatenating the vectors bLjk,bRjk,k=1,...,pj, respectively, as follows: 
BLj=[bLj1bLjpj],BRj=[bRj1bRjpj]
(84)
The corner points of the smooth folds in S0 defined in Eq. (38) provide a simple choice for the locations where the path γj(η) crosses the boundary rulings of the smooth folds. Thus, as in Sec. 3.3, the position vectors of such corner points are utilized to define the vectors bLjk,bRjk,k=1,...,pj (cf. Eq. (75)) 
[BLjBRj]=[12(|Cj|*+Cj)I312(|Cj|*Cj)I312(|Cj|*Cj)I312(|Cj|*+Cj)I3][P̂1P̂4]
(85)
After defining all the shape variables and geometric parameters associated with the smooth folds crossed by each path γj(η),j=1,...,NP, the folding map for points in the faces is obtained as the composition of transformations Hjk associated with the smooth folds crossed by these paths (cf. Eq. (63)) 
x=χ(X,t)where[x1]=(k=1pjHjk)[X1]
(86)
where Xspan(e1,e2) is the position vector of a point in P0jS0 and x3 is the position vector of such a point in PtjSt. Since the results of Theorem 3.1 and Eqs. (63) and (64) are also applicable for the smooth folds crossed by the paths γj(η), the matrices Hjk follow the same formulation presented in Eq. (64) 
Hjk=(T(bLjk)Q3(φ(mjk))Q1(ajkθjk)Q31(φ(mjk))T1(bLjk)T1((wjk0wjk)(e3×mjkmjk))T(bRjk)Q3(φ(mjk))Q1((1ajk)θjk)Q31(φ(mjk))T1(bRjk))
(87)
Determining the folding map for the smooth folds requires further steps since these surface subdomains undergo nonrigid deformations in addition to rotations and translations. The position vector of a point in Fti expressed in its associated fold coordinate system (with basis {ê1i,ê2i,ê3i} and origin at (1/2)(ci(1)+ci(1)), see Sec. 3.1) is denoted x̂3. Since the current configuration of a smooth fold Fti is known in its associated fold coordinate system for a given set of fold shape variables, the focus here is on the formulation of the map x=χ̂(x̂(t),t). To construct this map, a transformation represented by the matrix L̂i4×4 is first applied to the vector x̂. This transformation returns the position of the deformed smooth fold as joined to the reference configuration of the face adjacent to its boundary ruling Fi(1,ζ2) (denoted P0iF). The transformation L̂i is obtained by first aligning position of the smooth fold expressed in its associated fold coordinate system with the plane spanned by e1 and e2 (performed via the transformation Q1(âiθ̂i)T([0,ŵi/2,0])), and then aligning such a resulting position to P0iF (performed through the transformation T(p̂1i)Q3(φ(v̂i2v̂i1)). Therefore, L̂i is given as follows: 
L̂i=T(p̂1i)Q3(φ(v̂i2v̂i1))Q1(âiθ̂i)T([0,ŵi/2,0])
(88)
Subsequently, the folding map associated with the face P0iF is applied to determine the position of the smooth fold in the current configuration expressed in the global coordinate system. Thus, the folding map for points in the smooth folds is defined as follows: 
x=χ̂(x̂(t),t)where[x1]=(k=1piFH˘iFk)L̂i[x̂(t)1]
(89)
where x̂ is the position vector of a point in FtiSt expressed in the ith fold coordinate system with basis {ê1i,ê2i,ê3i} and x3 is the position vector of such a point expressed in the global coordinate system with basis {e1,e2,e3}. Note that the position vector of a point in F0iS0 is obtained as X=χ̂(x̂(0),0)span(e1,e2) 
X=χ̂(x̂(0),0)where[X1]=I4L̂i|t=0[x̂(0)1]=T(p̂1i)Q3(φ(v̂i2v̂i1))T([0,ŵi0/2,0])[x̂(0)1]
(90)

Numerical Implementation.

As stated at the beginning of Sec. 3, the continuous motion of the sheet is achieved by continuously altering the values of the fold shape variables (θ̂i,ŵi, and âi,i=1,...,NF) such that any attained configuration satisfies the constraints presented in Eq. (76) (i.e., the motion of the sheet is a continuous path in the constrained configuration space).

The simulation of the continuous motion of the sheet is executed by incrementally updating the values of the fold shape variables using input guess increments and then iteratively applying any required corrections such that the resulting set of fold shape variables satisfies the constraints of Eq. (76). The iterative procedure used to perform such corrections is outlined in this section. As stated in Sec. 3.3, the full set of fold shape variables θ̂i,ŵi, and âi,i=1,...,NF, is reduced to only the set of fold angles θ̂i,i=1,...,NF, by establishing relations ŵi(θ̂i,ŵi0) and âi(θ̂i,ŵi0) following certain assumptions on the extensibility and curvature field of the smooth folds (see Appendix  B for details).

Constraints.

Given a guess set of fold angles ordered in the vector lkθ̂NF where the subscript l denotes to incremental step number and the superscript k denotes to the correction iteration number, the matrices Rj(lkθ̂) and the vectors dj(lkθ̂),j=1,...,NI, are calculated. If the vector of fold angles lkθ̂ does not yield a configuration that satisfies Eq. (76), Rj(lkθ̂)I3 is not equal to the zero matrix in 3×3 and/or dj(lkθ̂) is not equal to 03.

Let R(lkθ̂)6NI+2NF with components j(lkθ̂) be the vector of residuals from constraints of Eq. (76) (6NI in total) and from constraints imposing the upper and lower fold angle bounds (2NF in total). Since Rj(lkθ̂) is an orthogonal matrix, only three of its scalar components are independent. Thus, the matrix-type constraint in Eq. (65) provides the following three scalar constraints that contribute to the residual vector [40]: 
6j5(lkθ̂)=12λR(R23j(lkθ̂))26j4(lkθ̂)=12λR(R31j(lkθ̂))26j3(lkθ̂)=12λR(R12j(lkθ̂))2
(91)
where j{1,...,NI} and λR is the weight for residuals from Eq. (65). This constant weight is an algorithmic parameter included in the residual vector to ensure that its components are scaled to a similar order of magnitude (i.e., the components of Rj(lkθ̂) are dimensionless while those of dj(lkθ̂) have units of length). The three components of the vector dj(lkθ̂),j=1,...,NI, which must be zero for the constraint in Eq. (66) to be satisfied, provide the following components to the residual vector: 
6j2(lkθ̂)=12λd(d1j(lkθ̂))26j1(lkθ̂)=12λd(d2j(lkθ̂))26j(lkθ̂)=12λd(d3j(lkθ̂))2
(92)

where j{1,...,NI} and λd is the weight for residuals from Eq. (66).

Fold Angle Bounds.

Additional components of R(lkθ̂) are included to ensure that the fold angles remain within prescribed upper and lower bounds during the continuous motion of the sheet. A penalty approach is used to implement the required fold angle bounds for each fold as in Ref. [40]. The lower bound for θ̂i is denoted θ̂iL[π,0] and the upper bound is denoted θ̂iU[0,π]. Conventional assignments for θ̂iL and θ̂iU are provided in Table 1.

The additional components of R(lkθ̂) required to enforce the lower bound of lkθ̂i consist of a penalty that is zero if lkθ̂iθ̂iL and increases proportionally to the square of the difference between lkθ̂i and θ̂iL when lkθ̂i<θ̂iL 
6NI+2i1(lkθ̂)=12λBmax(0,lkθ̂i+θ̂iL)2
(93)
where i{1,...,NF} and λB is the weight for residuals from fold angle bound constraints. Similarly, to enforce the upper bound of θ̂i, a penalty that is zero if lkθ̂iθ̂iU and increases proportionally to the square of the difference between lkθ̂i and θ̂iU when lkθ̂i>θ̂iU is included in the following components of R(lkθ̂): 
6NI+2i(lkθ̂)=12λBmax(0,lkθ̂iθ̂iU)2
(94)

where i{1,...,NF}.

Method of Solution.

Following the generalized Newton's method [40], the first-order expansion of R(θ̂) is used to determine the increment in θ̂ required to minimize the components of the constraint residual vector R(θ̂) 
R(θ̂+Δθ̂)=R(θ̂)+R(θ̂)θ̂Δθ̂+...=06NI+2NF
(95)
The following correction increment Δθ̂ is obtained from the previous first-order expansion: 
Δθ̂=(R(θ̂)θ̂)R(θ̂)
(96)

where (·) represents the Moore–Penrose pseudoinverse.

The lth set of guess fold angle increments is collected in the vector lINΔθ̂NF. First, the guess fold angle vector for the lth increment (l1θ̂) is calculated as follows: 
l1θ̂=l1θ̂+lΔθ̂
(97)
where the fold angles are first updated with the input fold angle increments lINΔθ̂ projected into the null space of the previous residual derivatives (see Ref. [40] for details) 
lΔθ̂=(INF(R(l1θ̂)θ̂)(R(l1θ̂)θ̂))lINΔθ̂
(98)
If R(l1θ̂)/(6NI+2NF)tol1, the fold angles are corrected iteratively as follows: 
lkΔθ̂=(R(lkθ̂)θ̂)R(lkθ̂)
(99)
 
lk+1θ̂=lkθ̂+lkΔθ̂
(100)

cf. Eq. (96). The correction process of Eqs. (99) and (100) is repeated until R(lk+1θ̂)/(6NI+2NF)<tol1 or lkΔθ̂/NF<tol2, where tol1 and tol2 are numerical tolerances. Table 2 summarizes the numerical procedure used for the kinematic simulation of sheets with smooth folds.

Implementation Examples

The numerical procedure used for the simulation of the deformation of sheets with smooth folds presented in Sec. 3.5 is implemented in matlab. The smooth folds Fti are visualized using the matlab three-dimensional shaded surface plot function surf while the faces Pti are visualized as filled three-dimensional polygons using fill3. The complete matlab code for the kinematic simulation of origami with smooth folds and various input examples are included in the “Supplemental Materials” tab for this paper on the ASME Digital Collection.

Smooth folds having G2 continuous joints with their adjacent faces are assumed for all the examples presented in this section. The formulation of the parametric curve ĉi(ζ1) that defines the cross-sectional shape of such smooth folds is provided in Eq. (34). Inextensible deformation and a symmetric parabolic curvature field are assumed for ĉi(ζ1) (see Appendix  B for details). Also, the range of fold angle values assumed for all the examples presented in this section is the interval [π,π]. It is noted that the aforementioned assumptions are taken for simplicity and do not present a limitation of the proposed model. Other assumptions on the extensibility and curvature field of the folds are also applicable as long as such assumptions do not violate the continuity conditions for ĉi(ζ1).

A sheet having eight smooth folds meeting at one interior fold intersection is shown in Fig. 13. The folds are enumerated in counterclockwise order. Various guess fold angle increments are considered ranging from simple to more complex. The folded configurations shown in Fig. 13(a) are obtained through the following guess fold angle increments:

 
lINΔθ̂=π50[10001000]l{1,...,50}
(101)
and the folded configurations shown in Fig. 13(b) are obtained through the following guess fold angle increments: 
lINΔθ̂=π50[01000100]l{1,...,50}
(102)
The two previous guess fold angle increments represent simple examples, and the fold angle correction procedure (refer to Table 2) converged prior to performing an initial correction iteration (i.e., R(l1θ̂)/NF<tol1l{1,...,50}). Alternatively, an example of a more complex folding deformation resulting from guess fold angle increments that required iterative corrections is shown in Fig. 13(c). For this example, the guess fold angle increments are as follows: 
lINΔθ̂=π50[11111111]l{1,...,50}
(103)

As shown in the fold angle versus increment plot in Fig. 13(c), the fold angles obtained from the simulation procedure differ from the simple addition of the guess fold angle increments. As observable in the configurations shown in Fig. 13(c), all the surface subdomains comprising the sheet remain joined through the deformation of the sheet. Thus, the simulation procedure presented in Sec. 3.5 successfully allows for the correction of the fold angles such that they satisfy the necessary constraints for a valid configuration presented in Eq. (76).

More complex examples of origami sheets having two interior fold intersections are shown in Fig. 14. The graphs of Fig. 15 show the vertices and fold centerlines associated with such sheets. Since the two interior fold intersections for these sheets share a common adjacent fold, their associated constraint equations are coupled. For all the three sheets, the guess fold angle increments are as follows (see Fig. 15 for the numbering of the folds):

 
lINΔθ̂=5π600[111111111111111]l{1,...,100}
(104)

The sheet with the baseline fold pattern shown in Fig. 14(a) exhibits a symmetric behavior, which is expected from the symmetry of the fold pattern and the guess fold angle increments. Sheets having fold patterns obtained by modifying the interior vertex coordinates and the boundary vertex coordinates of the baseline fold pattern are, respectively, shown in Figs. 14(b) and 14(c). It is observed that both sheets undergo dissimilar fold angle histories compared to the sheet having the baseline fold pattern as observed from both the folded configurations in Fig. 14 and the fold angle versus increment number plots in Fig. 16. These examples show the versatility of the present model that allows for simulation of sheets having arbitrary fold patterns and boundary shapes. This permits the potential application of the model in the design of origami-based structures and mechanisms constrained by realistic material and structural constraints.

The present model is also applicable for the simulation of sheets having arbitrary topology. To illustrate this, sheets having four, five, and six interior fold intersections are shown in Fig. 17.10 Various configurations are shown for these examples. As observed from these schematics, the present model captures well the behavior of the folded sheets during their full range of motion (fold angles vary from 0 to ±π for various folds in these sheets). The sheet having four interior fold intersections exhibits a desired response that the proposed model captures well, which is of configurations having all the faces located in parallel tangent planes but not overlapping one another. This may aid ongoing efforts to model and design rigid origami structures with thick sheets (see Refs. [12,8890]).

Summary and Conclusions

A model for the kinematic response of origami structures with smooth folds having nonzero sheet surface area and arbitrary order of geometric continuity was presented in this work. A brief review of an established model for rigid origami with conventional creased folds that is extended herein is provided in Sec. 2. The geometrical description of smooth folds was presented in Sec. 3.1, and the parametric representations of the fold cross-sectional shape for various orders of continuity were provided in Sec. 3.1.1. The fold pattern description (Sec. 3.2), the constraints on the sheet deformation for origami with smooth folds analogous to those for origami with creased folds (Sec. 3.3), and the map between reference and current configurations (Sec. 3.4) were also presented. The numerical implementation of the model allowing for simulation of the motion of sheets with arbitrary fold patterns was described in Sec. 3.5, and the implementation examples were provided in Sec. 4.

From the model formulation developed and the results presented in this work, a number of conclusions can be drawn. First, the proposed model successfully allows for the mathematical representation of origami with folds of nonzero sheet surface area and arbitrary order of geometric continuity (in terms of fold shape geometry, constraints on the fold shape variables, and mapping between reference and current configurations). The conventional rigid origami with creased folds of zeroth-order continuity represents a special case of this more general model and is captured as well. Furthermore, the arbitrary order of geometric continuity in the sheet considered in this work allows for the physical modeling of origami structures having significant thickness using plate or shell representations, which is a currently active research topic. Second, the present model and its associated numerical implementation allow for the simulation of sheets having arbitrary fold pattern geometry and topology and boundary shape. The careful consideration of constraints ensures that only meaningful valid configurations will be predicted. The proposed modeling and simulation framework readily allows for modifications to the fold pattern and folding sequence as demonstrated in the implementation examples provided in this work. These characteristics make such a framework useful in future fold pattern design and fold planning studies for origami structures and mechanisms having folds that cannot be accurately represented as creases.

Acknowledgment

This work was supported by the National Science Foundation and the Air Force Office of Scientific Research under Grant No. EFRI-1240483.

Appendix A: Derivation of the Constraint for Origami With Creased Folds

The constraint for origami with creased folds presented in Eq. (14) has been derived in the literature (see Refs. [36,37]) using an approach analogous to that used in Sec. 3.3 to derive the constraints for origami with smooth folds (Eqs. (65) and (66)). The purpose here is not to reproduce such a derivation from the literature but rather derive Eq. (14) as a special case of Eqs. (65) and (66) occurring when the smooth folds are degenerated to straight lines (corresponding to creased folds).

Since Eq. (14) is already accounted for in origami with smooth folds (cf. Eq. (65)), it is only needed to show that dj=03 (Eq. (66)) holds for origami with creased folds if the constraint presented in Eq. (14) is met. A segment of a general path γj(η) enclosing an interior fold intersection (corresponding to a single point for origami with creased folds) crossing each of its incident folds once in a counterclockwise order is shown in Fig. 18. Using Fig. 18 as a reference, the vectors w̃jk and l̃jk (see Eqs. (61) and (62)) are, respectively, given as follows:

 
w̃jk=[bjk00]
(A1)
 
l̃jk={[ljk+1cos(αjk)ljkbjkljk+1sin(αjk)0];k=1,...,nj1[lj1cos(αjk)ljkbjklj1sin(αjk)0];k=nj
(A2)

where ljk0,bjk, and ljk+bjk0.

Substituting Eqs. (A1) and (A2) into the expression for dj provided in Eq. (66) and utilizing the fact that R1(ϕ)w̃jk=w̃jk for the expression of w̃jk in Eq. (A1), the following is obtained: 
dj=k=1nj((l=1k1R1(θjl)R3(αjl))[bjk00])+(k=1nj1((l=1k1R1(θjl)R3(αjl))R1(θjk)[ljk+1cos(αjk)ljkbjkljk+1sin(αjk)0]))+((l=1nj1R1(θjl)R3(αjl))R1(θjnj)[lj1cos(αjnj)ljnjbjnjlj1sin(αjnj)0])=(k=1nj1((l=1k1R1(θjl)R3(αjl))R1(θjk)[ljk+1cos(αjk)ljkljk+1sin(αjk)0]))+((l=1nj1R1(θjl)R3(αjl))R1(θjnj)[lj1cos(αjnj)ljnjlj1sin(αjnj)0])
(A3)
The following equality is then utilized to simplify the previous expression: 
[ljk+1cos(αjk)ljk+1sin(αjk)0]=R3(αjk)R1(θjk+1)[ljk+100]k{1,...,nj1}
(A4)
Substituting Eq. (A4) into Eq. (A3), the following simplified expression is obtained: 
dj=R1(θj1)[lj100]+(l=1nj1R1(θjl)R3(αjl))R1(θjnj)[lj1cos(αjnj)lj1sin(αjnj)0]=[lj100]+(l=1nj1R1(θjl)R3(αjl))R1(θjnj)R3(αjnj)[lj100]=[lj100]+(l=1njR1(θjl)R3(αjl))[lj100]=[lj100]+Rj[lj100]
(A5)

The previous equation shows that Rj=I3dj=03 for origami with creased folds independently from the choice of the path γj(η). Therefore, dj=03 is a redundant constraint for origami with creased folds. This result shows that the constraints for origami with creased folds are a special case of the more general constraints for origami with smooth folds.

Appendix B: Determination of Fold Cross Section Shape Variables

As stated in Sec. 3.3, the full set of fold shape variables θ̂i,ŵi, and âi,i=1,...,NF, is reduced for simplification to only the set of fold angles θ̂i,i=1,...,NF, by determining relations ŵi=ŵi(θ̂i,ŵi0) and âi=âi(θ̂i,ŵi0). These relations are obtained by making assumptions on the extensibility and curvature field of the smooth folds. The process followed to determine such relations is outlined in this appendix.

Smooth folds having G2 continuous joints with their adjacent faces are considered. The formulation of the parametric curve ĉi(ζ1) describing the cross-sectional shape for these smooth folds is provided in Eq. (34). In addition to ŵi(θ̂i,ŵi0) and âi(θ̂i,ŵi0), the relations βL1i(θ̂i,ŵi0),βR1i(θ̂i,ŵi0),βL2i(θ̂i,ŵi0), and βR2i(θ̂i,ŵi0) are also determined here since they are required to fully define the shape of the smooth folds (see Eqs. (29) and (31)) even though these variables do not appear in the constraints presented in Sec. 3.3.

To determine the fold shape variables for any value of fold width ŵi0, the parametric curve ĉi(ζ1) is made nondimensional by ŵi0 as follows: 
c¯i(ζ1):=ĉi(ζ1)ŵi0
(B1)
This leads to the following nondimensional form of the parametric representation of ĉi(ζ1) presented in Eq. (34): 
c¯i(ζ1)=h50(ζ1)c¯L0i+h51(ζ1)c¯R0i+h52(ζ1)c¯L1i+h53(ζ1)c¯R1i+h54(ζ1)c¯L2i+h55(ζ1)c¯R2i
(B2)
where 
c¯L0i:=ĉL0iŵi0=[012ŵiŵi00],    c¯R0i:=ĉR0iŵi0=[012ŵiŵi00]
(B3)
 
c¯L1i:=ĉL1iŵi0=βL1iŵi0[0cos(âiθ̂i)sin(âiθ̂i)]c¯R1i:=ĉR1iŵi0=βR1iŵi0[0cos((1âi)θ̂i)sin((1âi)θ̂i)]
(B4)
 
c¯L2i:=ĉL2iŵi0=βL2iŵi0[0cos(âiθ̂i)sin(âiθ̂i)]c¯R2i:=ĉR2iŵi0=βR2iŵi0[0cos((1âi)θ̂i)sin((1âi)θ̂i)]
(B5)

Equations (B3)(B5) show that the nondimensional fold cross section c¯i(ζ1) is a function of the nondimensional variables θ̂i,ŵi/ŵi0,âi,βL1i/ŵi0,βR1i/ŵi0,βL2i/ŵi0, and βR2i/ŵi0.

The nondimensional arc length of c¯i(ζ1), denoted s¯itot, is determined as follows: 
s¯itot=sitotŵi0=11dc¯i(ζ1)dζ1dζ1
(B6)
and the nondimensional curvature of c¯i(ζ1), denoted κ¯i(ζ1), is given as follows: 
κ¯i(ζ1):=ŵi0κi(ζ1)=dc¯i(ζ1)dζ1×d2c¯i(ζ1)dζ12/dc¯i(ζ1)dζ13
(B7)
The nondimensional signed curvature of c¯i(ζ1), denoted κ¯is(ζ1), is given as 
κ¯is(ζ1)=ŵi0κis(ζ1)=κ¯i(ζ1)sgn((dc¯i(ζ1)dζ1×d2c¯i(ζ1)dζ12)·ê1i)
(B8)
The next step is to provide assumptions on the extensibility and curvature field of the fold cross section exhibited during folding. Inextensible deformation is assumed here. It should be noted that inextensibility of the folds is not a restriction of the present model but rather assumed for simplicity. Thus, the goal arc length of the fold cross section, sitotG, is constant 
sitotG=ŵi0
(B9)
A parabolic form for the goal signed curvature field κisG(ζ1) is assumed 
κisG(ζ1):=κimax(4si(ζ1)2(ŵi0)2+4si(ζ1)ŵi0)
(B10)

where si(ζ1) is defined in Eq. (23). Note that κisG(ζ1)=0 at si(ζ1)=0 and si(ζ1)=ŵi0.

The nondimensional form of Eq. (B9) yields the following nondimensional goal arc length for c¯i(ζ1), denoted s¯itotG: 
s¯itotG=sitotGŵi0=1
(B11)
The nondimensional form of the goal signed curvature field in Eq. (B10), denoted κ¯isG(ζ1), is the following: 
κ¯isG(ζ1)=ŵi0κisG(ζ1)=κ¯imax(4s¯i(ζ1)2+4s¯i(ζ1))
(B12)
where 
κ¯imax:=ŵi0κimax
(B13)
 
s¯i(ζ1)=si(ζ1)ŵi0=1ζ1dc¯i(ζ)dζdζ
(B14)
Once the assumptions on extensibility and curvature field are selected, the nondimensional variables θ̂i,ŵi/ŵi0,âi,βL1i/ŵi0,βR1i/ŵi0,βL2i/ŵi0, and βR2i/ŵi0 are fitted to satisfy such assumptions. For the symmetric smooth folds resulting from the selected assumptions, the variables to fit are simplified as follows: 
âi=12,βL1iŵi0=βR1iŵi0=:β1iŵi0,βL2iŵi0=βR2iŵi0=:β2iŵi0
(B15)
For a given value of κ¯imax, the error between the goal arc length (s¯itotG) and the arc length of c¯i(ζ1) (s¯itot), and the error between κ¯is and κ¯isG are minimized as follows: 
Findθ̂i[0,π],  ŵiŵi0(0,1]β1iŵi0(0,n],  β2iŵi0[0,n]thatminimizef(θ̂i,ŵiŵi0,β1iŵi0,β2iŵi0)
(B16)
where 
f(θ̂i,ŵiŵi0,β1iŵi0,β2iŵi0)=(s¯itotGs¯itot)2(s¯itotG)2+11(κ¯isG(ζ1)κ¯is(ζ1))2dζ111(κ¯isG(ζ1))2dζi
(B17)

It is assumed that the fold shape variables ŵi/ŵi0,β1i/ŵi0, and β2i/ŵi0 are equal for folds having the same fold angle absolute value |θ̂i|. Thus, the range for fold angles in the bounds of the optimization problem in Eq. (B16) only includes non-negative numbers.

The minimization problem in Eq. (B16) is repeated for a set of values of κ¯max to obtain values of θ̂i,ŵi,β1i,andβ2i associated with such a maximum nondimensional curvature, which satisfy the assumptions on extensibility and curvature field. The problem is solved for each value of κ¯max using the gradient-based optimization algorithm in matlabfmincon. The upper bound n for β1i/ŵi and β2i/ŵi is selected as 10, which is far from the actual values obtained for these parameters, and thus, it resulted in an inactive bound. Afterward, approximations for ŵi(θ̂i,ŵi0),β1i(θ̂i,ŵi0),β2i(θ̂i,ŵi0) and κimax(θ̂i,ŵi0) are interpolated from the data obtained by solving Eq. (B16) for a set of values for κ¯imax.

2

· denotes the two-norm, i.e., y=(y·y)1/2.

3

The matrix Kronecker product ⊗ is defined as [78]: YZ:{m×n,p×q}mp×nq, where [YZ]ijblockp×q=YijZ.

4

Self-intersection avoidance is an essential restriction in origami as stated in Definition 2.1. It remains an open problem to provide constraints on fold angles that would allow for three-dimensional folded configurations free of self-intersection [3,36]. This restriction is not considered in this work.

5

A ruled surface is formed by the union of straight lines, called the rulings or generators of the surface [79].

6

The sign function is defined as follows: sgn(y):={1;y<01;y>00;y=0.

7

As stated in Sec. 2.2, self-intersection avoidance is not considered in this work.

8

Alternative approaches for modeling the large rotations resulting from folding include quaternions-based [35] and geometric algebra-based approaches [87]. However, such approaches are not explored in this work.

9

A breve (˘) is used to distinguish the symbols related to the folding map from the symbols used in Sec. 3.3.

10

The example shown in Fig. 5 is also a simulation result obtained using the proposed model.

References

References
1.
Lang
,
R. J.
,
2007
, “
The Science of Origami
,”
Phys. World
,
20
(
2
), pp.
30
31
.
2.
Demaine
,
E. D.
,
2001
, “
Folding and Unfolding Linkages, Paper, and Polyhedra
,”
Discrete and Computational Geometry
,
Springer-Verlag
,
Berlin
, pp.
113
124
.
3.
Demaine
,
E. D.
, and
O’Rourke
,
J.
,
2007
,
Geometric Folding Algorithms
,
Cambridge University Press
,
Cambridge, UK
.
4.
Cromvik
,
C.
, and
Eriksson
,
K.
,
2009
, “
Airbag Folding Based on Origami Mathematics
,”
Fourth International Meeting of Origami Science, Mathematics, and Education
,
Origami 4
, Pasadena, CA, Sept. 8–10, pp.
129
139
.
5.
Pandey
,
S.
,
Gultepe
,
E.
, and
Gracias
,
D. H.
,
2013
, “
Origami Inspired Self-Assembly of Patterned and Reconfigurable Particles
,”
J. Visualized Exp.
,
72
, p.
e50022
.
6.
Morgan
,
J.
,
Magleby
,
S. P.
,
Lang
,
R. J.
, and
Howell
,
L. L.
,
2015
, “
A Preliminary Process for Origami-Adapted Design
,”
ASME
Paper No. DETC2015-47559.
7.
Gray
,
S.
,
Zeichner
,
N.
,
Kumar
,
V.
, and
Yim
,
M.
,
2011
, “
A Simulator for Origami-Inspired Self-Reconfigurable Robots
,”
Fifth International Meeting of Origami Science, Mathematics, and Education
,
Origami 5
, CRC Press, Taylor & Francis Group, Boca Raton, FL, pp.
323
333
.
8.
Hawkes
,
E.
,
An
,
B.
,
Benbernou
,
N. M.
,
Tanaka
,
H.
,
Kim
,
S.
,
Demaine
,
E. D.
,
Rus
,
D.
, and
Wood
,
R. J.
,
2010
, “
Programmable Matter by Folding
,”
Proc. Natl. Acad. Sci.
,
107
(
28
), pp.
12441
12445
.
9.
Lang
,
R. J.
,
2009
, “
Computational Origami: From Flapping Birds to Space Telescopes
,”
25th Annual Symposium on Computational Geometry
, Aarhus, Denmark, June 8–10, ACM, New York, pp.
159
162
.
10.
Wilson
,
L.
,
Pellegrino
,
S.
, and
Danner
,
R.
,
2013
, “
Origami Sunshield Concepts for Space Telescopes
,”
AIAA
Paper No. 2013-1594.
11.
Pohl
,
D.
, and
Wolpert
,
W. D.
,
2009
, “
Engineered Spacecraft Deployables Influenced by Nature
,”
Proc. SPIE
,
7424
, p.
742408
.
12.
Zirbel
,
S. A.
,
Lang
,
R. J.
,
Thomson
,
M. W.
,
Sigel
,
D. A.
,
Walkemeyer
,
P. E.
,
Trease
,
B. P.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2013
, “
Accommodating Thickness in Origami-Based Deployable Arrays
,”
ASME J. Mech. Des.
,
135
(
11
), p.
111005
.
13.
Cheng
,
Q.
,
Song
,
Z.
,
Ma
,
T.
,
Smith
,
B. B.
,
Tang
,
R.
,
Yu
,
H.
,
Jiang
,
H.
, and
Chan
,
C. K.
,
2013
, “
Folding Paper-Based Lithium-Ion Batteries for Higher Areal Energy Densities
,”
Nano Lett.
,
13
(
10
), pp.
4969
4974
.
14.
Nam
,
I.
,
Kim
,
G.-P.
,
Park
,
S.
,
Han
,
J. W.
, and
Yi
,
J.
,
2014
, “
All-Solid-State, Origami-Type Foldable Supercapacitor Chips With Integrated Series Circuit Analogues
,”
Energy Environ. Sci.
,
7
(
3
), pp.
1095
1102
.
15.
Song
,
Z.
,
Ma
,
T.
,
Tang
,
R.
,
Cheng
,
Q.
,
Wang
,
X.
,
Krishnaraju
,
D.
,
Panat
,
R.
,
Chan
,
C. K.
,
Yu
,
H.
, and
Jiang
,
H.
,
2014
, “
Origami Lithium-Ion Batteries
,”
Nat. Commun.
,
5
, p.
3140
.
16.
White
,
P. J.
,
Latscha
,
S.
,
Schlaefer
,
S.
, and
Yim
,
M.
,
2011
, “
Dielectric Elastomer Bender Actuator Applied to Modular Robotics
,”
IEEE/RSJ International Conference on Intelligent Robots and Systems
(
IROS
), San Francisco, CA, Sept. 25–30, pp.
408
413
.
17.
Lee
,
D.-Y.
,
Kim
,
J.-S.
,
Kim
,
S.-R.
,
Koh
,
J.-S.
, and
Cho
,
K.-J.
,
2013
, “
The Deformable Wheel Robot Using Magic-Ball Origami Structure
,”
ASME
Paper No. DETC2013-13016.
18.
Snyder
,
M. P.
,
Sanders
,
B.
,
Eastep
,
F. E.
, and
Frank
,
G. J.
,
2009
, “
Vibration and Flutter Characteristics of a Folding Wing
,”
J. Aircr.
,
46
(
3
), pp.
791
799
.
19.
Cheung
,
K. C.
,
Tachi
,
T.
,
Calisch
,
S.
, and
Miura
,
K.
,
2014
, “
Origami Interleaved Tube Cellular Materials
,”
Smart Mater. Struct.
,
23
(
9
), p.
094012
.
20.
Silverberg
,
J. L.
,
Evans
,
A. A.
,
McLeod
,
L.
,
Hayward
,
R. C.
,
Hull
,
T.
,
Santangelo
,
C. D.
, and
Cohen
,
I.
,
2014
, “
Using Origami Design Principles to Fold Reprogrammable Mechanical Metamaterials
,”
Science
,
345
(
6197
), pp.
647
650
.
21.
Schenk
,
M.
, and
Guest
,
S. D.
,
2013
, “
Geometry of Miura-Folded Metamaterials
,”
Proc. Natl. Acad. Sci.
,
110
(
9
), pp.
3276
3281
.
22.
Fuchi
,
K.
,
Diaz
,
A. R.
,
Rothwell
,
E. J.
,
Ouedraogo
,
R. O.
, and
Tang
,
J.
,
2012
, “
An Origami Tunable Metamaterial
,”
J. Appl. Phys.
,
111
(
8
), p.
084905
.
23.
Filipov
,
E. T.
,
Tachi
,
T.
, and
Paulino
,
G. H.
,
2015
, “
Origami Tubes Assembled Into Stiff, Yet Reconfigurable Structures and Metamaterials
,”
Proc. Natl. Acad. Sci.
,
112
(
40
), pp.
12321
12326
.
24.
Thrall
,
A. P.
, and
Quaglia
,
C. P.
,
2014
, “
Accordion Shelters: A Historical Review of Origami-Like Deployable Shelters Developed by the U.S. Military
,”
Eng. Struct.
,
59
, pp.
686
692
.
25.
Martínez-Martín
,
F. J.
, and
Thrall
,
A. P.
,
2014
, “
Honeycomb Core Sandwich Panels for Origami-Inspired Deployable Shelters: Multi-Objective Optimization for Minimum Weight and Maximum Energy Efficiency
,”
Eng. Struct.
,
69
, pp.
158
167
.
26.
Turner
,
N.
,
Goodwine
,
B.
, and
Sen
,
M.
, “
A Review of Origami Applications in Mechanical Engineering
,”
Proc. Inst. Mech. Eng., Part C
,
230
(
14
), pp.
2345
2362
.
27.
Peraza-Hernandez
,
E. A.
,
Hartl
,
D. J.
,
Malak
,
R. J.
, Jr.
, and
Lagoudas
,
D. C.
,
2014
, “
Origami-Inspired Active Structures: A Synthesis and Review
,”
Smart Mater. Struct.
,
23
(
9
), p.
094001
.
28.
Lebée
,
A.
,
2015
, “
From Folds to Structures, A Review
,”
Int. J. Space Struct.
,
30
(
2
), pp.
55
74
.
29.
Tachi
,
T.
,
2009
, “
Simulation of Rigid Origami
,”
Fourth International Meeting of Origami Science, Mathematics, and Education
,
Origami 4
, Pasadena, CA, Sept. 8–10, pp.
175
187
.
30.
Evans
,
T. A.
,
Lang
,
R. J.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2015
, “
Rigidly Foldable Origami Gadgets and Tessellations
,”
R. Soc. Open Sci.
,
2
(
9
), p.
150067
.
31.
Abel
,
Z.
,
Cantarella
,
J.
,
Demaine
,
E. D.
,
Eppstein
,
D.
,
Hull
,
T. C.
,
Ku
,
J. S.
,
Lang
,
R. J.
, and
Tachi
,
T.
,
2015
, “
Rigid Origami Vertices: Conditions and Forcing Sets
,”
e-print arXiv:1507.01644
.
32.
Tachi
,
T.
,
2010
, “
Geometric Considerations for the Design of Rigid Origami Structures
,”
International Association for Shell and Spatial Structures
(
IASS
) Symposium, Shanghai, China, Nov. 8–12, Vol.
12
, pp.
458
460
.
33.
Tachi
,
T.
,
2013
, “
Designing Freeform Origami Tessellations by Generalizing Resch’s Patterns
,”
ASME J. Mech. Des.
,
135
(
11
), p.
111006
.
34.
Tachi
,
T.
,
2013
, “
Freeform Origami Tessellations by Generalizing Resch’s Patterns
,”
ASME
Paper No. DETC2013-12326.
35.
Wu
,
W.
, and
You
,
Z.
,
2010
, “
Modelling Rigid Origami With Quaternions and Dual Quaternions
,”
Proc. R. Soc. London, Ser. A
,
466
(
2119
), pp.
2155
2174
.
36.
Belcastro
,
S.-M.
, and
Hull
,
T. C.
,
2002
, “
Modelling the Folding of Paper Into Three Dimensions Using Affine Transformations
,”
Linear Algebra Appl.
,
348
(
13
), pp.
273
282
.
37.
Belcastro
,
S.-M.
, and
Hull
,
T. C.
,
2002
, “
A Mathematical Model for Non-Flat Origami
,”
Third International Meeting of Origami Mathematics, Science, and Education
,
Origami 3
, pp.
39
51
.
38.
Tachi
,
T.
,
2007
, “
Rigid Origami Simulator
,”
TSG
, University of Tokyo College of Arts and Sciences, Komaba, Tokyo, Japan.
39.
Tachi
,
T.
,
2010
, “
Freeform Origami
,”
TSG
, University of Tokyo College of Arts and Sciences, Komaba, Tokyo, Japan.
40.
Tachi
,
T.
,
2010
, “
Freeform Variations of Origami
,”
J. Geom. Graphics
,
14
(
2
), pp.
203
215
.
41.
Schenk
,
M.
, and
Guest
,
S. D.
,
2011
, “
Origami Folding: A Structural Engineering Approach
,”
Fifth International Meeting of Origami Science, Mathematics, and Education
,
Origami 5
, pp.
291
304
.
42.
Demaine
,
E.
,
Demaine
,
M.
,
Koschitz
,
D.
, and
Tachi
,
T.
,
2011
, “
Curved Crease Folding: A Review on Art, Design and Mathematics
,”
IABSE-IASS Symposium: Taller, Longer, Lighter
(
IABSE-IASS2011
), London, Sept. 20–23, pp.
20
23
.
43.
Demaine
,
E. D.
,
Demaine
,
M. L.
,
Huffman
,
D. A.
,
Koschitz
,
D.
, and
Tachi
,
T.
,
2015
, “
Characterization of Curved Creases and Rulings: Design and Analysis of Lens Tessellations
,” Origami 6: I. Mathematics, pp.
209
230
.
44.
Dias
,
M. A.
,
Dudte
,
L. H.
,
Mahadevan
,
L.
, and
Santangelo
,
C. D.
,
2012
, “
Geometric Mechanics of Curved Crease Origami
,”
Phys. Rev. Lett.
,
109
(
11
), p.
114301
.
45.
Francis
,
K. C.
,
Rupert
,
L. T.
,
Lang
,
R. J.
,
Morgan
,
D. C.
,
Magleby
,
S. P.
, and
Howell
,
L. L.
,
2014
, “
From Crease Pattern to Product: Considerations to Engineering Origami-Adapted Designs
,”
ASME
Paper No. DETC2014-34031.
46.
Fuchi
,
K.
, and
Diaz
,
A. R.
,
2013
, “
Origami Design by Topology Optimization
,”
ASME J. Mech. Des.
,
135
(
11
), p.
111003
.
47.
Kergosien
,
Y. L.
,
Gotoda
,
H.
, and
Kunii
,
T. L.
,
1994
, “
Bending and Creasing Virtual Paper
,”
IEEE Comput. Graphics Appl.
,
14
(
1
), pp.
40
48
.
48.
Solomon
,
J.
,
Vouga
,
E.
,
Wardetzky
,
M.
, and
Grinspun
,
E.
,
2012
, “
Flexible Developable Surfaces
,”
Comput. Graphics Forum
,
31
(
5
), pp.
1567
1576
.
49.
Hwang
,
H.-D.
, and
Yoon
,
S.-H.
,
2015
, “
Constructing Developable Surfaces by Wrapping Cones and Cylinders
,”
Comput.-Aided Des.
,
58
, pp.
230
235
.
50.
Schreck
,
C.
,
Rohmer
,
D.
,
Hahmann
,
S.
,
Cani
,
M.-P.
,
Jin
,
S.
,
Wang
,
C. C.
, and
Bloch
,
J.-F.
,
2015
, “
Nonsmooth Developable Geometry for Interactively Animating Paper Crumpling
,”
ACM Trans. Graphics
,
35
(
1
), p.
10
.
51.
Zhu
,
L.
,
Igarashi
,
T.
, and
Mitani
,
J.
,
2013
, “
Soft Folding
,”
Comput. Graphics Forum
,
32
(
7
), pp.
167
176
.
52.
Fuchi
,
K.
,
Ware
,
T. H.
,
Buskohl
,
P. R.
,
Reich
,
G. W.
,
Vaia
,
R. A.
,
White
,
T. J.
, and
Joo
,
J. J.
,
2015
, “
Topology Optimization for the Design of Folding Liquid Crystal Elastomer Actuators
,”
Soft Matter
,
11
(
37
), pp.
7288
7295
.
53.
Fuchi
,
K.
,
Buskohl
,
P. R.
,
Ware
,
T.
,
Vaia
,
R. A.
,
White
,
T. J.
,
Reich
,
G. W.
, and
Joo
,
J. J.
,
2014
, “
Inverse Design of LCN Films for Origami Applications Using Topology Optimization
,”
ASME
Paper No. SMASIS2014-7497.
54.
Peraza Hernandez
,
E. A.
,
Hu
,
S.
,
Kung
,
H. W.
,
Hartl
,
D.
, and
Akleman
,
E.
,
2013
, “
Towards Building Smart Self-Folding Structures
,”
Comput. Graphics
,
37
(
6
), pp.
730
742
.
55.
Peraza-Hernandez
,
E. A.
,
Hartl
,
D. J.
, and
Malak
,
R. J.
, Jr.
,
2013
, “
Design and Numerical Analysis of an SMA Mesh-Based Self-Folding Sheet
,”
Smart Mater. Struct.
,
22
(
9
), p.
094008
.
56.
Peraza-Hernandez
,
E.
,
Hartl
,
D.
,
Galvan
,
E.
, and
Malak
,
R.
,
2013
, “
Design and Optimization of a Shape Memory Alloy-Based Self-Folding Sheet
,”
ASME J. Mech. Des.
,
135
(
11
), p.
111007
.
57.
Peraza Hernandez
,
E.
,
Hartl
,
D.
,
Malak
,
R.
, and
Lagoudas
,
D.
,
2015
, “
Analysis and Optimization of a Shape Memory Alloy-Based Self-Folding Sheet Considering Material Uncertainties
,”
ASME
Paper No. SMASIS2015-9001.
58.
Ahmed
,
S.
,
Ounaies
,
Z.
, and
Frecker
,
M.
,
2014
, “
Investigating the Performance and Properties of Dielectric Elastomer Actuators as a Potential Means to Actuate Origami Structures
,”
Smart Mater. Struct.
,
23
(
9
), p.
094003
.
59.
McGough
,
K.
,
Ahmed
,
S.
,
Frecker
,
M.
, and
Ounaies
,
Z.
,
2014
, “
Finite Element Analysis and Validation of Dielectric Elastomer Actuators Used for Active Origami
,”
Smart Mater. Struct.
,
23
(
9
), p.
094002
.
60.
Liu
,
Y.
,
Boyles
,
J. K.
,
Genzer
,
J.
, and
Dickey
,
M. D.
,
2012
, “
Self-Folding of Polymer Sheets Using Local Light Absorption
,”
Soft Matter
,
8
(
6
), pp.
1764
1769
.
61.
Mailen
,
R. W.
,
Liu
,
Y.
,
Dickey
,
M. D.
,
Zikry
,
M.
, and
Genzer
,
J.
,
2015
, “
Modelling of Shape Memory Polymer Sheets That Self-Fold in Response to Localized Heating
,”
Soft Matter
,
11
(
39
), pp.
7827
7834
.
62.
Davis
,
D.
,
Mailen
,
R.
,
Genzer
,
J.
, and
Dickey
,
M. D.
,
2015
, “
Self-Folding of Polymer Sheets Using Microwaves and Graphene Ink
,”
RSC Adv.
,
5
(
108
), pp.
89254
89261
.
63.
Ionov
,
L.
,
2013
, “
Nature-Inspired Stimuli-Responsive Self-Folding Materials
,”
Intelligent Stimuli-Responsive Materials: From Well-Defined Nanostructures to Applications
,
Wiley
,
New York
, pp.
1
16
.
64.
Demaine
,
E. D.
, and
Demaine
,
M. L.
,
2002
, “
Recent Results in Computational Origami
,”
Third International Meeting of Origami Mathematics, Science, and Education
,
Origami 3
, pp.
3
16
.
65.
Fuchi
,
K.
,
Buskohl
,
P. R.
,
Bazzan
,
G.
,
Durstock
,
M. F.
,
Reich
,
G. W.
,
Vaia
,
R. A.
, and
Joo
,
J. J.
,
2015
, “
Origami Actuator Design and Networking Through Crease Topology Optimization
,”
ASME J. Mech. Des.
,
137
(
9
), p.
091401
.
66.
Lang
,
R. J.
,
1996
, “
A Computational Algorithm for Origami Design
,”
12th Annual Symposium on Computational Geometry
, Philadelphia, PA, May 24–26, ACM, New York, pp.
98
105
.
67.
Tachi
,
T.
,
2010
, “
Origamizing Polyhedral Surfaces
,”
IEEE Trans. Visualization Comput. Graphics
,
16
(
2
), pp.
298
311
.
68.
Peraza Hernandez
,
E.
,
Hartl
,
D.
,
Malak
,
R.
,
Akleman
,
E.
,
Gonen
,
O.
, and
Kung
,
H.
,
2016
, “
Design Tools for Patterned Self-Folding Reconfigurable Structures Based on Programmable Active Laminates
,”
ASME J. Mech. Rob.
,
8
(
3
), p.
031015
.
69.
Peraza-Hernandez
,
E.
,
Frei
,
K.
,
Hartl
,
D.
, and
Lagoudas
,
D.
,
2014
, “
Folding Patterns and Shape Optimization Using SMA-Based Self-Folding Laminates
,”
Proc. SPIE
,
9057
, p.
90571G
.
70.
Zhou
,
X.
,
Wang
,
H.
, and
You
,
Z.
,
2015
, “
Design of Three-Dimensional Origami Structures Based on a Vertex Approach
,”
Proc. R. Soc. London, Ser. A
,
471
(
2181
), p.
20150407
.
71.
Saito
,
K.
,
Tsukahara
,
A.
, and
Okabe
,
Y.
,
2016
, “
Designing of Self-Deploying Origami Structures Using Geometrically Misaligned Crease Patterns
,”
Proc. R. Soc. London, Ser. A
,
472
(
2185
), p.
20150235
.
72.
Song
,
G.
, and
Amato
,
N. M.
,
2004
, “
A Motion-Planning Approach to Folding: From Paper Craft to Protein Folding
,”
IEEE Trans. Rob. Autom.
,
20
(
1
), pp.
60
71
.
73.
Amato
,
N. M.
,
Dill
,
K. A.
, and
Song
,
G.
,
2003
, “
Using Motion Planning to Map Protein Folding Landscapes and Analyze Folding Kinetics of Known Native Structures
,”
J. Comput. Biol.
,
10
(
3–4
), pp.
239
255
.
74.
Xi
,
Z.
, and
Lien
,
J.-M.
,
2014
, “
Folding Rigid Origami With Closure Constraints
,”
ASME
Paper No. DETC2014-35556.
75.
Xi
,
Z.
, and
Lien
,
J.-M.
,
2015
, “
Folding and Unfolding Origami Tessellation by Reusing Folding Path
,”
IEEE International Conference on Robotics and Automation
(
ICRA
), Seattle, WA, May 26–30, pp.
4155
4160
.
76.
Xi
,
Z.
, and
Lien
,
J.-M.
,
2015
, “
Plan Folding Motion for Rigid Self-Folding Machine Via Discrete Domain Sampling
,”
IEEE International Conference on Robotics and Automation
(
ICRA
), Seattle, WA, May 26–30, pp.
2938
2943
.
77.
Kuribayashi
,
K.
,
Tsuchiya
,
K.
,
You
,
Z.
,
Tomus
,
D.
,
Umemoto
,
M.
,
Ito
,
T.
, and
Sasaki
,
M.
,
2006
, “
Self-Deployable Origami Stent Grafts as a Biomedical Application of Ni-Rich TiNi Shape Memory Alloy Foil
,”
Mater. Sci. Eng.: A
,
419
(
12
), pp.
131
137
.
78.
Horn
,
R. A.
, and
Johnson
,
C. R.
,
1991
,
Topics in Matrix Analysis
,
Cambridge University Press
,
Cambridge, UK
.
79.
Pressley
,
A. N.
,
2010
,
Elementary Differential Geometry
,
Springer-Verlag, London
.
80.
Calladine
,
C. R.
,
1989
,
Theory of Shell Structures
,
Cambridge University Press
,
Cambridge, UK
.
81.
Akleman
,
E.
, and
Chen
,
J.
,
2006
, “
Insight for Practical Subdivision Modeling With Discrete Gauss–Bonnet Theorem
,”
International Conference on Geometric Modeling and Processing
, Pittsburgh, PA, July 26–28, Springer-Verlag, Berlin, pp. 287–298.
82.
Sullivan
,
J. M.
,
2008
, “
Curvatures of Smooth and Discrete Surfaces
,”
Discrete Differential Geometry
,
Birkhauser
,
Basel, Switzerland
, pp.
175
188
.
83.
Tachi
,
T.
,
2015
, “
Rigid Folding of Periodic Origami Tessellations
,”
Origami 6: I. Mathematics
, American Mathematical Society, Providence, RI, pp.
97
108
.
84.
Barsky
,
B. A.
, and
DeRose
,
T. D.
,
1984
, “Geometric Continuity of Parametric Curves,”
Computer Science Division, University of California
,
Berkeley, CA
,
Technical Report No. UCB/CSD 84/205
.
85.
Barsky
,
B. A.
,
Bartels
,
R. H.
, and
Beatty
,
J. C.
,
1987
,
An Introduction to Splines for Use in Computer Graphics and Geometric Modeling
,
M. Kaufmann Publishers
,
Los Altos, CA
.
86.
Cheney
,
W.
, and
Kincaid
,
D.
,
1996
,
Numerical Analysis. Mathematics of Scientific Computing
,
Brooks & Cole Publishing Company
,
Pacific Grove, CA
.
87.
McRobie
,
F.
, and
Lasenby
,
J.
,
2000
, “
The Kinematics of Large Rotations Using Clifford Algebra
,”
IUTAM-IASS
Symposium on Deployable Structures: Theory and Applications, Springer Science+Business Media Dordrecht, The Netherlands, pp.
271
280
.
88.
Ku
,
J. S.
, and
Demaine
,
E. D.
,
2015
, “
Folding Flat Crease Patterns With Thick Materials
,”
ASME
Paper No. DETC2015-48039.
89.
Chen
,
Y.
,
Peng
,
R.
, and
You
,
Z.
,
2015
, “
Origami of Thick Panels
,”
Science
,
349
(
6246
), pp.
396
400
.
90.
Tachi
,
T.
,
2011
, “
Rigid-Foldable Thick Origami
,”
Fifth International Meeting of Origami Science, Mathematics, and Education
,
Origami 5
, CRC Press, Boca Raton, FL, pp.
253
264
.

Supplementary data