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 [9–12], electronic components with improved properties [13–15], robotic components [16,17], foldable wings [18], cellular materials [19], metamaterials [20–23], and shelters [24,25], among others [26–28].

*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,32–34]. 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 (*G*^{0}) 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. [42–44]); 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 [47–51]. 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 [47–51] 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 [54–57], dielectric elastomer [58,59], and optically responsive polymeric self-folding sheets [60–62], 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,64–71] 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. [72–76] 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 $ei\u2208\mathbb{R}3$, *i* = 1, 2, 3, with $e3:=e1\xd7e2$ 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 $P0i\u2282S0,i=1,...,NP$, where $NP$ is the number of faces in the sheet (i.e., $S0=\u222ai=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 $Pti\u2282St,i=1,...,NP$, i.e., $St=\u222ai=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 configuration*$St$*has 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 *G*^{0} 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*$\theta \u0302i(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 $\theta \u0302i(t),\u2009i=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 $\theta \u0302i$ on *t* is left implicit to simplify the notation.

### 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 in*$S0$*. Each vertex has an associated position vector denoted*$vj\u2208span(e1,e2)$.

*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 $C\u0302\u2208{\u22121,0,1}NF\xd7(NI+NB)$ with elements $C\u0302ij$ is introduced

Let *n _{j}*, $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 $mjk\u2208span(e1,e2),\u2009j=1,...,NI,\u2009k=1,...,nj$, be the vector along the length of the

*k*th fold line incident to the

*j*th interior vertex that emanates from such a vertex (see Fig. 4).

*The angle between the vector*$y\u2208span(e1,e2)$

*and*$e1$

*,*

*which starts at*$e1$

*and is measured in the counterclockwise direction, is denoted*$\phi (y):span(e1,e2)\u2192[0,2\pi )$

*and is determined as follows*:

^{2}

*j*th interior vertex and are defined as follows:

^{3}

where $In$ denotes the $\mathbb{R}n\xd7n$ identity matrix.

### Constraints.

^{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

*K*, is considered [54,80,81]. It is given as $2\pi $ 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

_{j}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\pi $ 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

*α*as defined in $S0$. As a consequence, developability is assured for any valid current configuration $St$.

_{jk}Nontrivial and important are the constraints on the fold angles that define the constrained configuration space. The set of constraints for $\theta \u0302i,\u2009i=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].

*k*th fold incident to the

*j*th interior vertex is denoted

*θ*. The vectors $\theta j\u2208\mathbb{R}nj,\u2009j=1,...,NI$, are formed by collecting the fold angles $\theta jk,\u2009k=1,...,nj,$ as follows:

_{jk}where the elements of the matrix $Cj$ are defined in Eq. (6), and $|\xb7|*:\mathbb{R}m\xd7n\u2192\mathbb{R}\u22650m\xd7n$ denotes the elementwise absolute value of a matrix where $[|Y|*]ij=|Yij|$.

The fold angles associated with the folds incident to each interior vertex $\theta jk,\u2009k=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:

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*$\theta j1=0$. - (ii)
*If the jth interior vertex has two incident creased folds, it allows for a valid configuration if (ii,1)*$\alpha j1\u2260\pi ,\theta j1=\theta j2=0$*or (ii,2)*$\alpha j1=\pi ,\theta j1=\theta j2$. - (iii)
*If the jth interior vertex has three incident creased folds, it allows for a valid configuration if (iii,1)*$\alpha j1\u2260\pi ,\alpha j2\u2260\pi ,\alpha j3\u2260\pi ,\theta j1=\theta j2=\theta j3=0$*, or (iii,2)*$\alpha j1=\pi ,\theta j1=\theta j2,\theta j3=0$*, or (iii,3)*$\alpha j2=\pi ,\theta j2=\theta j3,\theta j1=0$*, or (iii,4)*$\alpha j3=\pi ,\theta j3=\theta j1,\theta j2=0$.

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

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 $ei\u2208\mathbb{R}3,\u2009i=1,2,3$, with $e3:=e1\xd7e2$ 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,\u2009i=1,...,NP,\u2009F0i,\u2009i=1,...,NF$, and $I0i,\u2009i=1,...,NI$, respectively. Therefore, $S0=(\u222ai=1NPP0i)\u222a(\u222ai=1NFF0i)\u222a(\u222ai=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,\u2009i=1,...,NP,\u2009Fti,\u2009i=1,...,NF$, and $Iti,\u2009i=1,...,NI$, respectively. Thus, $St=(\u222ai=1NPPti)\u222a(\u222ai=1NFFti)\u222a(\u222ai=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 configuration*$St$*has 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.

*Smooth folds*:

*The smooth folds are ruled surfaces*

^{5}of the following form:*where*$Fti(\zeta 1,\zeta 2)\u2208\mathbb{R}3$*is a parameterization of*$Fti$*. Without loss of generality, the parameters ζ _{1} and ζ_{2} are contained in the intervals*$[\u22121,1]$

*and*$[0,1]$

*, respectively.*

An example of a smooth fold is shown in Fig. 6. In Definition 3.2, $hti\u2208\mathbb{R}3$ provides the direction of the rulings comprising $Fti$ while $cti(\zeta 1)\u2208\mathbb{R}3$ is the parametric ruled surface directrix curve that defines the cross section of $Fti$. The curve parameterized by $cti(\zeta 1)$ is contained in a plane orthogonal to $hti$ as stated in Eq. (16). It is assumed that $\u2225hti\u2225$ 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(\zeta 1)$. To simplify the notation, the dependence of $Fti(\zeta 1,\zeta 2),\u2009cti(\zeta 1)$, and $hti$ on *t* is left implicit for the remainder of the paper and these vectors are denoted as $Fi(\zeta 1,\zeta 2),\u2009ci(\zeta 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(\zeta 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(\u22121,\zeta 2)$ and $Fi(1,\zeta 2)$, refer to Fig. 6). The remaining boundaries of $Fti$ (i.e., $Fi(\zeta 1,0)$ and $Fi(\zeta 1,1)$) are either at the boundary of $St$ or joined to a fold intersection. A parameterization for the fold intersections ($Iti,\u2009i=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*$\theta \u0302i(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 $\theta \u0302\u2208\mathbb{R}NF$ was constructed by collecting the fold angles $\theta \u0302i,\u2009i=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(\zeta 1)$ are also provided.

A schematic of a smooth fold cross section showing $w\u0302i$ is provided in Fig. 7. The vector $w\u0302\u2208\mathbb{R}NF$ is constructed by collecting the variables $w\u0302i$ as follows:

*fold width*$w\u0302i0$ is the value of $w\u0302i$ at

*t*= 0. The vector $w\u03020\u2208\mathbb{R}NF$ is constructed as follows:

The fold-attached orthonormal vectors $e\u0302ji\u2208\mathbb{R}3,\u2009i=1,...,NF$, *j* = 1, 2, 3, with $e\u03023i:=e\u03021i\xd7e\u03022i$ form the bases ${e\u03021i,e\u03022i,e\u03023i}$ that define the local *fold coordinate system* of each smooth fold $Fti$. The origin of this coordinate system is located at $(1/2)(ci(\u22121)+ci(1))$. The director vector $hi$ is aligned to $e\u03021i$ (i.e., $hi\xb7e\u03021i=\u2225hi\u2225$) and $(ci(1)\u2212ci(\u22121))\u2208span(e\u03022i)$.

^{6}

*G*

^{0}continuous joints with the faces adjacent to $Fti$ require the following conditions on $c\u0302i(\zeta 1)$ at $\zeta 1=\xb11$ (refer to Fig. 7):

*G*

^{1}continuous joints with the planar faces adjacent to $Fti$ in addition to

*G*

^{0}continuity [84,85]. The following values of $ti(\zeta 1)$ at $\zeta 1=\xb11$ are then required for

*G*

^{1}continuity in addition to those conditions of Eq. (27):

*G*

^{1}continuity at the joints with the planar faces adjacent to $Fti$ :

where $\beta L1i,\u2009\beta R1i\u2208\mathbb{R}>0$.

*G*

^{2}continuity [85] in addition to

*G*

^{1}continuity. This requires the curvature of $c\u0302i(\zeta 1)$ at $\zeta 1=\xb11$ 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):

*G*

^{2}continuity at the joints with the planar faces adjacent to $Fti$ :

#### Fold Shape Examples.

*G*

^{1}continuity at the boundary rulings of $Fti,\u2009c\u0302i(\zeta 1)$ is expressed as follows:

*G*

^{2}continuity at the boundary rulings of $Fti,\u2009c\u0302i(\zeta 1)$ is expressed as follows:

### 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 in*$S0$*. Each vertex has an associated position vector denoted*$vj\u2208span(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 $\u2202S0$ 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 $\u2202S0$ 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).

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

As in Sec. 2.1, let *n _{j}*, $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 $mjk\u2208span(e1,e2),\u2009j=1,...,NI,\u2009k=1,...,nj$, be the vector along the length of the

*k*th fold centerline incident to the

*j*th interior vertex that emanates from this vertex. The vectors $mjk,\u2009k=1,...,nj$, have a counterclockwise ordering, i.e., $\phi (mj1)<\phi (mj1)<\cdots <\phi (mjnj)\u2200\u2009j\u2208{1,...,NI}$, where $\phi (\xb7)$ is defined in Eq. (4).

The mapping from the vertex position vectors $vj,\u2009j=1,...,NI+NB$, to the vectors $mjk,\u2009k=1,...,nj$, is provided in Eq. (7). The angles between adjacent fold centerlines intersecting at a common interior vertex $\alpha jk,\u2009j=1,...,NI,\u2009k=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 valid^{7} (according to Definition 3.1).

*α*, $k=1,...,nj$, satisfy the following constraint since they are defined in $S0$ :

_{jk}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 $c\u0302i(\zeta 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 $\theta \u0302i$, the distance $w\u0302i$ between the end-points of the cross section parametric curve $c\u0302i(\zeta 1)$, and $a\u0302i,\u2009i=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.

*θ*,

_{jk}*w*, and

_{jk}*a*are those associated with the

_{jk}*k*th smooth fold adjacent to $I0j$. The vectors $\theta j\u2208\mathbb{R}nj,\u2009j=1,...,NI$, constructed by collecting the fold angles

*θ*, $k=1,...,nj$, are defined in Eq. (10). The vectors $wj,\u2009aj\u2208\mathbb{R}nj,\u2009j=1,...,NI$, are constructed by collecting the variables

_{jk}*w*and

_{jk}*a*, $k=1,...,nj$, and are defined as follows:

_{jk}Note that if the vector $mjk$ is coincident with and has the *same* orientation as the *i*th fold centerline, then $ajk=a\u0302i$. Conversely, if the vector $mjk$ is coincident with and has the *opposite* orientation as the *i*th fold centerline, then $ajk=1\u2212a\u0302i$. This is obtained from the mapping provided in Eq. (44).

Let $\gamma j(\eta ):[0,1]\u2192S0$ be any simple closed path (i.e., $\gamma j(0)=\gamma j(1)$) enclosing $I0j$ and crossing each smooth fold adjacent to $I0j$ once. An example of a path $\gamma j(\eta )$ is shown in Fig. 10. The point having position $\gamma j(0)=\gamma 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 $\gamma j(\eta )$ 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 $\gamma j(\eta )$ crosses each boundary ruling of the smooth folds with associated vector $mjk$ are denoted $bLjk\u2208span(e1,e2)$ (point where $\gamma j(\eta )$ enters the smooth fold) and $bRjk\u2208span(e1,e2)$ (point where $\gamma j(\eta )$ exits the smooth fold). This is shown schematically in Fig. 10(a).

*b*

_{i}*The transformation matrix*$Ljk$

*describing the deformation associated with the folding of the kth smooth fold crossed by the path*$\gamma j(\eta )$

*is determined as follows*:

^{8}

*where the vectors*$gjk\u2208span(e1,e2),\u2009k=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:*

*Proof*. Referring to Figs. 11(c) and 11(d), the transformation associated with the folding of the *k*th smooth fold crossed by $\gamma j(\eta )$ can be discretized by two consecutive rotation transformations. The first transformation corresponds to a rotation by $(1\u2212ajk)\theta jk$ about an axis aligned to $mjk$ and crossing a point with position vector $bRjk\u2212gjk$. Using Eq. (50), the transformation matrix associated with this rotation is the following:

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 $\u2212(wjk0\u2212wjk)$. To define the position of the axes of rotation for the transformations composing $Ljk$, the vector $\u2211l=1k\u22121(wjl0\u2212wjl)[e3\xd7(mjl/\u2225mjl\u2225)]$ must be subtracted from $bLjk$ while $\u2211l=1k(wjl0\u2212wjl)[e3\xd7(mjl/\u2225mjl\u2225)]$ must be subtracted from $bRjk$. The recursive definition of the vectors $gjk,\u2009k=0,...,nj$, in Eq. (52) allows for a simplified form of these vector subtractions.□

*fixed*in space (not translating or rotating) for the derivation of constraints. Let $X\u2208span(e1,e2)$ be the position vector of a point in a face adjacent to $I0j$ in the reference configuration and let $x\u2208\mathbb{R}3$ be the position vector of such a point in a current configuration. The map $X\u21a6x$ is constructed as the composition of transformations $Ljk$ associated with the folds crossed by the segment of path $\gamma j(\eta )$ that connects $\gamma j(0)$ to the face containing the point with initial position

**X**

and *n _{y}* is the number of smooth folds crossed by the segment of the path $\gamma j(\eta )$ that connects $\gamma 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}*a*, and

_{jk}*w*, $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:

_{jk}*For the initially closed strip of faces and smooth folds joined to*$I0j$

*to remain closed with each face undergoing a rigid deformation, two constraints must hold. These are the rotation constraint*

*and the translation constraint*

*Proof*. The deformation map of a point in the face containing the point with position $\gamma j(0)$ must be the identity transformation since such a face is assumed fixed (i.e., if

*n*=

_{y}*n*in Eq. (63) is considered, then $x=X$). This requires the following:

_{j}cf. Eq. (66). The 21 and the 22 blocks of the right side of Eq. (67) are equal to $03\u22a4$ and 1, respectively.□

Corollary 3.1.

- (i)
*If*$I0j$*has a single adjacent smooth fold, it allows for a valid configuration if*$\theta j1=0$. - (ii)
*If*$I0j$*has two adjacent smooth folds, it allows for a valid configuration if (ii,1)*$\alpha j1\u2260\pi ,\theta j1=\theta j2=0$*or (ii,2)*$\alpha j1=\pi ,\theta j1=\theta j2$. - (iii)
*If*$I0j$*has three adjacent smooth folds, it allows for a valid configuration if (iii,1)*$\alpha j1\u2260\pi ,\alpha j2\u2260\pi ,\alpha j3\u2260\pi ,\theta j1=\theta j2=\theta j3=0$*, or (iii,2)*$\alpha j1=\pi ,\theta j1=\theta j2,\theta j3=0$*, or (iii,3)*$\alpha j2=\pi ,\theta j2=\theta j3,\theta j1=0$*, or (iii,4)*$\alpha j3=\pi ,\theta j3=\theta j1,\theta 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.

where the vectors $P\u0302i$ are defined in Eq. (40). Note that the mapping provided in Eq. (75) results in $bLjk=p\u03021i,\u2009bRjk=p\u03024i$ if $mjk$ is coincident with and has the same orientation as the *i*th smooth fold centerline and $bLjk=p\u03024i,\u2009bRjk=p\u03021i$ if $mjk$ is coincident with and has the opposite orientation as the *i*th smooth fold centerline.

cf. Eq. (15).

In general, the variables $\theta \u0302i,\u2009w\u0302i$, and $a\u0302i$ 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 $\theta \u0302i,\u2009w\u0302i$, and $a\u0302i$ 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 $c\u0302i(\zeta 1)$ are taken such that the overall deformation of a smooth fold becomes a function of the fold angle $\theta \u0302i$ and the fold width $w\u0302i0$ (i.e., $w\u0302i=w\u0302i(\theta \u0302i,w\u0302i0),\u2009a\u0302i=a\u0302i(\theta \u0302i,w\u0302i0)$ for each fold). Such a process involves nondimensionalization of the parametric curve $c\u0302i(\zeta 1)$ and is presented in detail in Appendix B for the case of smooth folds exhibiting *G*^{2} 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 $\gamma \u2323j(\eta ):[0,1]\u2192S0,\u2009j=1,...,NP$, be paths connecting the fixed face to $P0j,\u2009j=1,...,NP$, respectively (see Fig. 12 for an example). The paths $\gamma \u2323j(\eta ),\u2009j=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 $\gamma \u2323j(\eta )$ between $\gamma \u2323j(0)$ (located at the fixed face) and $\gamma \u2323j(1)$ (located at $P0j$). Each path $\gamma \u2323j(\eta )$ crosses a number of *p _{j}* smooth folds.

*positively*if it enters such a fold from the ruling $Fi(\u22121,\zeta 2)$ and exits at the ruling $Fi(1,\zeta 2)$. Conversely, the path $\gamma \u2323j(\eta )$ crosses a smooth fold

*negatively*if it crosses such a fold in the opposite direction. The matrices $C\u2323j\u2208{\u22121,0,1}pj\xd7NF,\u2009j=1,...,NP$, with elements $C\u2323kij$ are used for the identification and ordering of the folds crossed by $\gamma \u2323j(\eta )$

*k*th smooth fold crossed by $\gamma \u2323j(\eta )$. Such a vector has the same orientation as the fold centerline if $\gamma \u2323j(\eta )$ crosses the fold positively and opposite orientation if $\gamma \u2323j(\eta )$ crosses the fold negatively. The vector $M\u2323j\u2208\mathbb{R}3pj$ is constructed by concatenating the vectors $m\u2323jk,\u2009k=1,...,pj$, as follows:

*i*th fold coordinate system with basis ${e\u03021i,e\u03022i,e\u03023i}$ and $x\u2208\mathbb{R}3$ 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 $F0i\u2282S0$ is obtained as $X=\chi \u0302(x\u0302(0),0)\u2208span(e1,e2)$

### 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 ($\theta \u0302i,\u2009w\u0302i$, and $a\u0302i,\u2009i=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 $\theta \u0302i,\u2009w\u0302i$, and $a\u0302i,\u2009i=1,...,NF$, is reduced to only the set of fold angles $\theta \u0302i,\u2009i=1,...,NF$, by establishing relations $w\u0302i(\theta \u0302i,w\u0302i0)$ and $a\u0302i(\theta \u0302i,w\u0302i0)$ 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 $\u2009lk\theta \u0302\u2208\mathbb{R}NF$ where the subscript *l* denotes to incremental step number and the superscript *k* denotes to the correction iteration number, the matrices $Rj(lk\theta \u0302)$ and the vectors $dj(lk\theta \u0302),\u2009j=1,...,NI$, are calculated. If the vector of fold angles $\u2009lk\theta \u0302$ does not yield a configuration that satisfies Eq. (76), $Rj(lk\theta \u0302)\u2212I3$ is not equal to the zero matrix in $\mathbb{R}3\xd73$ and/or $dj(lk\theta \u0302)$ is not equal to $03$.

*λ*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\theta \u0302)$ are dimensionless while those of $dj(lk\theta \u0302)$ have units of length). The three components of the vector $dj(lk\theta \u0302),\u2009j=1,...,NI$, which must be zero for the constraint in Eq. (66) to be satisfied, provide the following components to the residual vector:

_{R}where $j\u2208{1,...,NI}$ and *λ _{d}* is the weight for residuals from Eq. (66).

#### Fold Angle Bounds.

Additional components of $R(lk\theta \u0302)$ 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 $\theta \u0302i$ is denoted $\theta \u0302iL\u2208[\u2212\pi ,0]$ and the upper bound is denoted $\theta \u0302iU\u2208[0,\pi ]$. Conventional assignments for $\theta \u0302iL$ and $\theta \u0302iU$ are provided in Table 1.

*λ*is the weight for residuals from fold angle bound constraints. Similarly, to enforce the upper bound of $\theta \u0302i$, a penalty that is zero if $\u2009lk\theta \u0302i\u2264\theta \u0302iU$ and increases proportionally to the square of the difference between $\u2009lk\theta \u0302i$ and $\theta \u0302iU$ when $\u2009lk\theta \u0302i>\theta \u0302iU$ is included in the following components of $R(lk\theta \u0302)$:

_{B}where $i\u2208{1,...,NF}$.

#### Method of Solution.

where $(\xb7)\u2020$ represents the Moore–Penrose pseudoinverse.

*l*th set of guess fold angle increments is collected in the vector $\u2009lIN\Delta \theta \u0302\u2208\mathbb{R}NF$. First, the guess fold angle vector for the

*l*th increment ($\u2009l1\theta \u0302$) is calculated as follows:

cf. Eq. (96). The correction process of Eqs. (99) and (100) is repeated until $\Vert R(\u2009lk+1\theta \u0302)\Vert /(6NI+2NF)<tol1$ or $\Vert lk\Delta \theta \u0302\Vert /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 *G*^{2} continuous joints with their adjacent faces are assumed for all the examples presented in this section. The formulation of the parametric curve $c\u0302i(\zeta 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 $c\u0302i(\zeta 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 $[\u2212\pi ,\pi ]$. 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 $c\u0302i(\zeta 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:

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):

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 $\xb1\pi $ 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,88–90]).

## 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 $\gamma j(\eta )$ 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\u0303jk$ and $l\u0303jk$ (see Eqs. (61) and (62)) are, respectively, given as follows:

where $ljk\u2208\mathbb{R}\u22650,\u2009bjk\u2208\mathbb{R}$, and $ljk+bjk\u22650$.

The previous equation shows that $Rj=I3\u21d2dj=03$ for origami with creased folds independently from the choice of the path $\gamma j(\eta )$. 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 $\theta \u0302i,\u2009w\u0302i$, and $a\u0302i,\u2009i=1,...,NF$, is reduced for simplification to only the set of fold angles $\theta \u0302i,\u2009i=1,...,NF$, by determining relations $w\u0302i=w\u0302i(\theta \u0302i,w\u0302i0)$ and $a\u0302i=a\u0302i(\theta \u0302i,w\u0302i0)$. 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 *G*^{2} continuous joints with their adjacent faces are considered. The formulation of the parametric curve $c\u0302i(\zeta 1)$ describing the cross-sectional shape for these smooth folds is provided in Eq. (34). In addition to $w\u0302i(\theta \u0302i,w\u0302i0)$ and $a\u0302i(\theta \u0302i,w\u0302i0)$, the relations $\beta L1i(\theta \u0302i,w\u0302i0),\u2009\beta R1i(\theta \u0302i,w\u0302i0),\u2009\beta L2i(\theta \u0302i,w\u0302i0)$, and $\beta R2i(\theta \u0302i,w\u0302i0)$ 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.

Equations (B3)–(B5) show that the nondimensional fold cross section $c\xafi(\zeta 1)$ is a function of the nondimensional variables $\theta \u0302i,\u2009w\u0302i/w\u0302i0,\u2009a\u0302i,\u2009\beta L1i/w\u0302i0,\u2009\beta R1i/w\u0302i0,\u2009\beta L2i/w\u0302i0$, and $\beta R2i/w\u0302i0$.

where $si(\zeta 1)$ is defined in Eq. (23). Note that $\kappa isG(\zeta 1)=0$ at $si(\zeta 1)=0$ and $si(\zeta 1)=w\u0302i0$.

It is assumed that the fold shape variables $w\u0302i/w\u0302i0,\u2009\beta 1i/w\u0302i0$, and $\beta 2i/w\u0302i0$ are equal for folds having the same fold angle absolute value $|\theta \u0302i|$. 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 $\kappa \xafmax$ to obtain values of $\theta \u0302i,\u2009w\u0302i,\u2009\beta 1i,\u2009and\beta 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 $\kappa \xafmax$ using the gradient-based optimization algorithm in matlabfmincon. The upper bound *n* for $\beta 1i/w\u0302i$ and $\beta 2i/w\u0302i$ 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 $w\u0302i(\theta \u0302i,w\u0302i0),\beta 1i(\theta \u0302i,w\u0302i0),\u2009\beta 2i(\theta \u0302i,w\u0302i0)$ and $\kappa imax(\theta \u0302i,w\u0302i0)$ are interpolated from the data obtained by solving Eq. (B16) for a set of values for $\kappa \xafimax$.

$\Vert \xb7\Vert $ denotes the two-norm, i.e., $\Vert y\Vert =(y\xb7y)1/2$.

The matrix Kronecker product ⊗ is defined as [78]: $Y\u2297Z:{\mathbb{R}m\xd7n,\mathbb{R}p\xd7q}\u2192\mathbb{R}mp\xd7nq$, where $[Y\u2297Z]ij\u2009block\u2208\mathbb{R}p\xd7q=YijZ.$

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.

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

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

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

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

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