## Abstract

With specific fold patterns, a 2D flat origami can be converted into a complex 3D structure under an external driving force. Origami inspires the engineering design of many self-assembled and re-configurable devices. This work aims to apply the level set-based topology optimization to the generative design of origami structures. The origami mechanism is simulated using thin shell models where the deformation on the surface and the deformation in the normal direction can be simplified and well captured. Moreover, the fold pattern is implicitly represented by the boundaries of the level set function. The folding topology is optimized by minimizing a new multiobjective function that balances kinematic performance with structural stiffness and geometric requirements. Besides regular straight folds, our proposed model can mimic crease patterns with curved folds. With the folding curves implicitly represented, the curvature flow is utilized to control the complexity of the folds generated. The performance of the proposed method is demonstrated by the computer generation and physical validation of two thin shell origami designs.

## 1 Introduction

### 1.1 Origami-Inspired Design: State of the Art.

Demonstrating the ability to transform from 2D to 3D, origami-inspired deployable structures have been applied in engineering [1], such as packaging, deployable solar panels, self-deployable mechanisms, and metamaterial design [2–5]. Mathematicians examined origami from a geometric perspective to comprehend the relationship between the fold pattern and the folded status [6]. Hull [7,8] studied the relation between the folding angles and the crease patterns and delivered a theorem for the local flat foldability. However, the theorem is insufficient for global foldability. Robert Lang applied the circle packing technique and proposed the software treemaker to compute the crease patterns from stick figures [9].

As for the analysis of the distortion of origami structures, various models have been developed. One well-adopted approach is the so-called rigid origami, which assumes that all the deformation occurs in the folding lines, and the facets between folding lines are rigid [8,10,11]. Tachi [10,12] developed simulation software for rigid origami. The configuration of the model, or the deformed shape, is explicitly represented by crease angles, and the trajectory is calculated by projecting angle motion into the constrained space. This pioneering work provided some freedom to the user to revise the folding lines and avoid local self-intersections. However, the global self-intersection still exists. Similarly, Wei et al. [13] studied the periodic pleated origami using geometric mechanics. To address the problem of nonnegligible fold thickness or with maximum curvature at the folds [14], Lagoudas and coworkers [15,16] presented a smooth folds model for rigid origami simulation, which can simulate the realistic folds of nonzero surface area and contain higher-order geometric continuity.

However, the rigid origami model cannot completely express the actual behavior of an origami. In reality, the facet area undergoes nonnegligible deformation, such as bending, stretching, and shearing. To address the problem, Schenk and Guest [11] proposed a pin-joint framework to track the motion of a partially folded sheet with a low computational cost. Their model also links the pure kinematics to the stiffness matrix approach. Liu and Paulino [17,18] proposed a nonlinear formulation to simulate large-displacement origami structures. They simplified the bar-and-hinge model [19] as pin-jointed bar networks with virtual rotational springs for origami simulations and thus reduced the degree-of-freedom (DOF) and improved the computational efficiency. The model can globally capture the essential features of origami, including folding, bending, and multistability. Based on the bar-and-hinge model, Ghassaei et al. [20] introduced a high-speed origami simulator using a graphics processing unit (GPU). This work is promising regarding the computational speed and interactivity that offer real-time simulation and interaction. However, compromises have to be considered between geometric accuracy and simulation speed, and instability will occur on facets with high aspect ratio triangles. Another engineering approach is using shell or membrane elements in the finite element analysis to simulate the deployment of origami [21,22]. Cai et al. [22] conducted the nonlinear finite element analysis using the variable Poisson’s ratio model to revise the stress distribution of membrane elements and successfully tracked the deployment of the famous Miura-Ori pattern.

Recently, in the field of computer graphics and architectures, the studies of origami with curved folds are emerging [23–25], since its potential in kinetic systems for architectural application [26] and shape-programmable structures [27,28]. In designing the compliant mechanism, Nelson et al. [29] addressed the importance of curved origami in designing developable structures. Their subsequent works [30–32] presented an energy method based on the normalized coordinate equations to identify the natural configurations of general curve folds. Moreover, in engineering, one of the significant applications of flexible electronics—the wearable soft electronics, compliant origami with curved creases—is more desirable because of easier conformal to the human body than conventional rigid origami.

Conventionally, origami design is based on the analysis and alteration of several well-known crease patterns. Fuchi et al. [33,34] did some pioneering work on topology optimization of origami structures. They proposed an optimization framework to achieve an optimal origami design by adding or removing folding lines to or from a reference crease pattern. Later, they enriched the work with a nonlinear truss model to model the sequenced motion and the large deformation of an origami [35,36]. In addition, a recent work proposed by Paulino and coworkers [37] offered a novel perspective of generating origami design. They applied a shape grammar formalism coupled with an interpreter to construct and modify an origami tessellation. These origami design frameworks are mainly based on rigid origami, composed of rigid patches with straight creases working as rotational joints.

### 1.2 Topology Optimization.

Topology optimization aims to find the best design geometry for optimal performance under certain constraints. The density-based approach, such as the homogenization method [38–40] and the solid isotropic material with penalization method [41,42], is the most popular way to perform topology optimization on the surface. However, the design obtained by the density-based optimization method can contain checkerboard patterns or gray elements. In contrast to the density-based method, the level set method can provide a clear boundary design. Moreover, since the level set functions are defined in the space with one higher dimension, the higher-order geometric information, such as curvatures and normal vectors, is embedded naturally in the geometric model. It allows the level set method for an exclusive capability of dealing with topological changes [43]. The level set-based topology optimization (TO) approach has been considered a powerful tool in generating innovative designs ever since the shape sensitivity analysis was cast into the framework [44–47].

This work presents a systematic solution to simulate and conceive origami-inspired compliant structures. We integrate the shell model with the level set topology optimization framework. The origami is modeled utilizing the finite element analysis with shell elements. Thus, the in-plane membrane, out-of-plane bending, and shear deformation can be well captured. Moreover, by incorporating the level set-based topology optimization method, the crease patterns are represented implicitly by the boundaries of the level set functions. The optimization problem is formulated as a multi-objective problem that aims to achieve an optimal distribution of the folding lines to fulfill the prescribed requirements. The topology of the crease pattern is optimized by solving the partial differential equation (PDE), the so-called Hamilton–Jacobi equation. Rather than adopting straight folding lines, our optimization method can form curved folds.

The remainder of this article is organized as follows: Sec. 2 introduces the background on the level set-based topology optimization framework and the shell model for simulating origami structures. In Sec. 3, we formulate the problem of optimizing origami structures and provide the sensitivity analysis. Next, an example of an origami gripper and an origami twisting mechanism is presented in Sec. 4. Section 5 first discusses the effects of initial designs and then shows the performance of the objective function through a numerical example. Finally, in Sec. 6, the conclusions are drawn, and the future work is also briefly discussed.

## 2 Method Overview

In this study, we propose a solution to systematic modeling and optimization of the origami structures to meet the mechanical requirements. The method is based on the level set topology optimization framework and a modified shell model. The flowchart is shown in Fig. 1.

### 2.1 Level Set Representation for Origami Structures.

*φ*are defined in

*R*

^{2}or

*R*

^{3}as an implicit function on one higher dimension [44]. The boundaries of each level set function represent the mountain and valley folds, respectively, and the areas apart from the boundary are the facets as shown in Eq. (1).

*D*is a bounded area represents the design domain, and

*D*⊂ R.

**x**is a point inside the design domain.

*k*represents the different level set functions. In our case, the number of

*k*is 2.

*M*and

*V*represent the mountain folds and valley folds, respectively, as shown in Fig. 2. The crease pattern’s topology is optimized by solving the Hamilton–Jacobi equation as Eq. (2), which is defined by differentiating the level set function with respect to time

*t*[44].

*φ*to represent a general level set function unless otherwise stated.

The presented Hamilton–Jacobi equation can be solved using the upwind spatial derivative to update the design iteratively [45]. A velocity extension step is needed to expand the normal velocity from the boundaries to the whole computational domain. In our case, the natural extension method is applied since the strain field is defined on the entire design domain [48].

Instead of assuming straight folding lines, this mode can form a general fold pattern with regular straight folds as long as curved ones. With the level set method, higher-order geometric information, such as the curvature, can be used to control the length and straightness of the generated crease patterns.

### 2.2 The Shell Model for Simulating Origami Structures.

*t*is the local thickness and

*E*represents Young’s modules. Moreover, we modify the local offset plane in the area near the crease with a sign relevant to the crease type, as shown in Fig. 3. The fold locations are determined by the boundary of the level set functions calculated by the following smoothed delta function [45].

*ɛ*controls the size of the transition area, which is chosen as two times the mesh size [47,43]. Then, we normalize

*δ*(

*φ*) to $\delta ^(\varphi )$, where the maximum value is equal to 1, and couple the $\delta ^(\varphi )$ on the shell model. The scheme of the level set representation of a modified shell model is shown in Fig. 3. It is worth noting that the smoothed delta function can provide a gradual transaction in the thickness direction. The middle surface location

**r**of the shell model can be represented as follows [50,51]:

**r**

*is the meshed surface,*

_{R}*ξ*

_{1}and

*ξ*

_{2}are the two curvilinear coordinates in the middle surface, and

*n*is the normal direction on the shell. The

*ζ*

_{0}is the offset in the thickness direction, which is related to $\delta ^(\varphi )$ as shown in Eq. (6).

*α*ranges from 1 to 2. Thus, with the offset, an additional term is added onto the strain and acts as a small perturbation to the model. We testify the modified shell model with different crease patterns with both straight and curved folds, and their deformed shapes are shown in Fig. 4. It is worth noting that the physics properties of the intersections between mountain and valley folding lines have already been implied. The intersections share the same thickness and Young’s modulus as folding lines, but their offsets in the thickness direction (

*ζ*

_{0}) have been canceled out. Thus, the overlaps are essentially areas with weak material but without the in-thickness offsets.

In the finite element analysis of the continuum shell, we adopt the classic Mindlin–Reissner plate shell theory that takes into account the through-thickness shear deformation [52–54]. The background information regarding the shell model and the derivation of the linear elastic strain energy density is presented in Appendix A. In terms of the numerical implementation, the SC6R element is applied, which consists of six nodes, and each node has six DOFs [53–55].

## 3 Problem Formulation and Shape Sensitivity Analysis

The topology optimization of an origami mechanism is equivalent to finding the optimum crease pattern to meet the engineering design requirements. In this work, we solve the optimization problem by minimizing a multi-objective function, which comprises three sections:

the engineering target (

*J*_{1}) to meet the kinematic requirements;the characteristic requirements for the design to be an origami structure (

*J*_{2});the perimeter constraint for achieving concise designs.

*U*is the space kinematically admissible displacement [55], Ω is the region occupied with the linear elastic material.

_{ad}*I*

_{1},

*I*

_{2}, and

*α*are constant weighting factors for kinematic target, characteristic objective, and perimeter constraint, respectively. The internal virtual energy

*a*(

**,**

*u***) and the exterior virtual energy**

*v**l*(

**) are defined as follows:**

*v**D*is the design domain; Γ = ∂Ω is the boundary of the design;

*a*(

**,**

*u***) is a symmetric bilinear function in terms of the displacement**

*v**u,*and the test function

*v*, that is,

*a*(

**,**

*u***) =**

*v**a*(

**,**

*v***); and**

*u**l*is a linear function in terms of the body force

*f*and the traction force

*g*on the Neumann boundary conditions Γ

*, as illustrated in Fig. 5. We treat each of the objectives as a subproblem and can make sensitivity analysis for each subproblem and collate the results. It is known that the sensitivity analysis result of the perimeter constraint is the mean curvature*

_{N}*κ*and can be calculated by [56]:

Therefore, we will only show the detailed steps for sensitivity analysis of the kinematic (*J*_{1}) and the characteristic requirements (*J*_{2}).

### 3.1 Shape Sensitivity Analysis of the Kinematic Objective.

*u*and a given target displacement

*u*

_{0}of a selected zone

*w*(

*x*) in the design domain;

*w*(

*x*) is zero excepted on the selected area.

*t*is formulated as follows:

### 3.2 Shape Sensitivity Analysis of the Characteristic Objective.

Since an origami structure shows the characteristic of the deformation concentrated on the folding lines, the facets are mainly flat, and the system undergoes large out-of-plane deformation [19]. Thus, we remount the objective functions by optimizing a structure, which, under a load scenario, the deformation of the system will be in favor of the following requirements:

The deformation occurs mainly in the crease area.

*J*_{2I}.The out-of-plane deformation dominates.

*J*_{2II}.

*J*

_{2I}and

*J*

_{2II}can be represented as following equations:

*a*(

**,**

*u***) is the total strain energy,**

*u**a*

_{Γ}(

**,**

*u***) is the total strain energy on the creases, and**

*u**a*(

_{b}**,**

*u***) is the total bending energy.**

*u**J*

_{2}as follows:

*J*

_{2}can be represented as follows:

*D*

_{Ω}() = ∂()/∂Ω(). In the following sections, we formulate the objective into two subproblems and derive the sensitivity for

*D*

_{Ω}

*J*

_{2I}and

*D*

_{Ω}

*J*

_{2II}separately.

#### 3.2.1 Minimizing the Deformation in the Facet Area.

*J*

_{2I}, which minimizes the difference of total strain energy

*a*(

**,**

*u***) and the total strain energy on the creases**

*u**a*

_{Γ}(

**,**

*u***), as shown in Eq. (22):**

*u***is the adjoint variable and is achieved by solving Eq. (27) using the finite element method (FEM).**

*v*#### 3.2.2 Maximizing the Out-of-Plane Deformation.

*J*

_{2II}can be formulated to maximize the bending energy on the crease patterns. With the observation that in the FEM model, the bending energy is concentrated on the area near the crease, we further simplify the objective function as maximum the bending energy on the design domain. As such, the objective function can be formulated as follows:

**, we solve the adjoint equation:**

*v*## 4 Numerical Implementation

### 4.1 An Origami Gripper Design.

*I*

_{1},

*I*

_{2}, and

*α*are constant weights of the kinematic target, characteristic objective, and perimeter constraint, respectively. To reduce the difficulty of adjusting the weighting factors, a linear relationship is set as

*I*

_{2}= 1−

*I*

_{2}. Thus, the velocity field can be represented as follows:

The design domain is a 1-by-1 square, with a thickness of 0.003. The boundary conditions are shown in Fig. 6(a). In the *x*-direction, the input displacements are set as 0.1, as shown in Fig. 6(a). The center point is fixed. The output areas *w*(*x*) are 0.1-by-0.05 rectangles in the middle of the top and bottom sides. The target displacement *u*_{0} is set to (0,0,1). Due to the symmetry, the upper left quarter of the domain is discretized with a rectangular 100 × 100 mesh. Young’s modulus is 10^{6} Pa on the facet and 5 × 10^{5} Pa on the folding lines.

The evolution of the 2D design is shown in Fig. 6(b), and the corresponding deformed status is shown in Fig. 6(c). Initially, the design is a flat thin shell consisting of circular valley folds. It can be observed that the design evolves from a group of uniformly distributed circles to a concise crease pattern. The optimization history plot is shown in Fig. 7. The plots fluctuate during the first 50 iterations, where the topology of the design changes significantly; then, after 200 iterations, the perimeter constraint dominates, so that it helps smooth the folding lines and the objective curve changes to stability. Finally, both objective plots converge around the 300th iteration.

### 4.2 An Origami Twister Design.

This example shows a twister design. The design domain is set to be a 1-by-0.5 rectangle, and its thickness is 0.003. The objective function consists of the mechanic requirements, the characteristic part, and the mean curvature control. The target of the design is set to achieve a twist motion on the tip areas of the rectangle. The boundary conditions are shown in Fig. 9(a). On the two side edges, the prescribed displacement in the *x*-direction is set as 0.1. The center area of the design domain is fixed. The output areas *w*(*x*) are 0.2-by-0.05 rectangles in the middle of the top and bottom sides. The target displacement on the top window is set to (0.1,0,0.1) and (−0.1,0,0.1) on the bottom. The asymmetric boundary condition is applied on the horizontal centerline of the rectangle. The upper domain is selected as the computational domain and discretized with a rectangular 100 × 50 mesh. Young’s modulus is 10^{6} Pa on the facet and 5 × 10^{5} Pa on the folding lines.

The evolution of the 2D design and their deformed status are shown in Fig. 9. The initial design is two groups of circular folding lines evenly distributed in the design domain. We can observe that the final design is bending outward, which is satisfied with the target. Figure 10 plots the displacement distribution in the *x*-direction, and the final design is shown on the lower side of the figure.

## 5 Discussion on the Initial Designs and the Characteristic Objective

### 5.1 The Effects of the Initial Design.

To test the method's robustness, we solve the problem again using different initial designs but with the same boundary conditions. Figure 11 shows the evolution of the design for two other models of different initial fold patterns. In contrast to the example shown in Fig. 6, the two designs introduce both the valley and the mountain folds as design variables. In detail, the top group has the same total number of initial folds as Fig. 6, and the bottom group is initially made up of more folding lines. It is worth noting that different initial designs can converge to the same optimum structure with the same imposed boundary conditions. In other words, the proposed method is robust in the sense that it shows little dependency on the initial designs, as shown in Fig. 11.

### 5.2 The Performance of the Characteristic Objective.

*J*

_{2}. The objective function only contains the

*J*

_{2}with the linear elastic constraint, as shown in the equations:

The design domain is a unit square. Figure 12 shows the applied boundary conditions (Fig. 12(a)), the initial design (Fig. 12(b)), and the achieved result (Fig. 12(c)). The center of the square is fixed, and two prescribed displacements *δ _{d}* in the

*x*-direction are assigned on the two sides to push the model move toward the center; the displacements in other directions are set as zero. The initial design consists of mountain and valley circular folding lines. With the effect of the objective, the design converges from circles to rhombuses. Since the perimeter constraints are not applied, the final result is not smooth and contains zigzags.

Figure 13 shows the evaluation of the performance of the given example. Figure 13(a) shows the convergence plot of each part of the objective function. The solid line plots the total strain energy on the facet area, and the hidden line shows the one over total bending energy in each iteration number. The two plots decrease in the first 20 iterations and then reach convergence and vibrate within a range. Figure 13(b) shows the strain energy density on the final design. We can observe that the highest number of energy concentrates on the folding lines, which means the deformation mainly occurs on the folding lines.

## 6 Conclusions

The purpose of this study is to provide a systematic way of modeling and optimizing an origami-like structure to meet engineering requirements, for instance, transforming force or displacement from one port to another. The origami is physically modeled by a shell model; thus, the in-plane stretching, out-of-plane bending, and shear deformation can be well captured.

Moreover, the folding lines are represented implicitly by the zero contours of the level set functions. Instead of assuming straight folding lines, our model can form curved folding lines, and the straightness of creases can be controlled by the perimeter constraint. An origami gripper and a twisting mechanism are optimized as demonstrations. The proposed method shows little dependence on the initial crease patterns as long as the boundary conditions are the same.

However, some aspects need to be improved. First, in the current setup of the problem, the design is based on linear elasticity and small deformation. In the future, we will examine the relatively large nonlinear deformation to capture the origami transformation with greater precision. Second, the current formulation of the problem relies mainly on maximizing or minimizing specific energies on the structure. We do not explicitly control the overlapping between mountain and valley folds. The evolution of the design is driven by the velocity field derived from the physics equation. However, one particular type of overlap, or vertex, should be discussed. Mathematicians [6] have shown that the number of vertices, mountain, and valley folding lines determines an origami’s foldability. The geometrical concerns, such as foldability requirements, need to be coupled into the problem formulation in the future work. In addition, the origami structure is optimized as a mechanism triggered by input loads in this work. The behaviors such as self-locking, face overlapping, and folding sequences have to be taken into account.

In terms of applications, the work can be expanded into the design of programmable folding surfaces by considering the manufacture of origami with active materials, such as ferromagnetic particles [57], thermal active polymers, and dielectric elastomers. Under a specific stimulus, it automatically switches between the original state and a target form. Moreover, with the conformal topology optimization framework with the extended level set method [57–62], we can optimize conformal origami structures on the surface with desired properties and kinematic performance.

## Acknowledgment

This project is supported by the National Science Foundation (CMMI-1762287), the Ford University Research Program (URP) (Award No. 2017-9198R), and the start-up fund from the State University of New York at Stony Brook.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

The data and information that support the findings of this article are freely available at: http://me.eng.stonybrook.edu/~chen/

### Appendix: The Strain Energy Density of a Linear Elastic Shell Structure

*ξ*

_{1}and

*ξ*

_{2}be the two curvilinear coordinates in the middle surface and

*ξ*

_{3}be the coordinate in the thickness direction varying within [−1 1] [50,51]. Thus, any point on the shell structure can be represented as

*(*

**x***ξ*

_{1}

*, ξ*

_{2}

*, ξ*

_{3}) [63,64].

*represents the reference coordinates. The energy bilinear form for linear elasticity is shown as follows:*

_{m}**stands for the displacement and**

*u***is the test function.**

*v**ξ*

_{1}

*,ξ*

_{2}coordinates. Thus, Eq. (42) can be simplified as follows [66]:

*ɛ*

_{ij}

^{1}(

**) and**

*u**ɛ*

_{ij}

^{2}(

**) are the membrane-shear stain and the bending strain, respectively.**

*u**ξ*

_{3}∈ [−1, 1], we can rewrite the energy bilinear form as follows [55,66]:

*ξ*

_{3}∈ [−1, 1], and thus, Eq. (46) can be simplified as follows [55]:

*a*(

_{b}**,**

*u***) to denote the bending strain energy in the contents mentioned earlier.**

*u*