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.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 , 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 . 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 .
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.  studied the periodic pleated origami using geometric mechanics. To address the problem of nonnegligible fold thickness or with maximum curvature at the folds , 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  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  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.  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.  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  and shape-programmable structures [27,28]. In designing the compliant mechanism, Nelson et al.  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  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 . 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.
The presented Hamilton–Jacobi equation can be solved using the upwind spatial derivative to update the design iteratively . 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 .
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.
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 (J1) to meet the kinematic requirements;
the characteristic requirements for the design to be an origami structure (J2);
the perimeter constraint for achieving concise designs.
Therefore, we will only show the detailed steps for sensitivity analysis of the kinematic (J1) and the characteristic requirements (J2).
3.1 Shape Sensitivity Analysis of the Kinematic Objective.
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 . 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. J2I.
The out-of-plane deformation dominates. J2II.
3.2.1 Minimizing the Deformation in the Facet Area.
3.2.2 Maximizing the Out-of-Plane Deformation.
4 Numerical Implementation
4.1 An Origami Gripper Design.
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 u0 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 106 Pa on the facet and 5 × 105 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 106 Pa on the facet and 5 × 105 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.
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.
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  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 , 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.
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/