Variational Level Set Method for Topology Optimization of Origami Fold Patterns

The flowchart

2.1 Level Set Representation for Origami Structures.

In this work, a level set representation for topology optimization of an origami structure is proposed. Basically, we divide the design domain into two sections: the crease (folding line) and the facet area. A convex (concave) crease is called the mountain, or valley, fold line. Thus, two level set functions are introduced. Similar to the conventional level set method, the level set functions φ are defined in R² or R³ 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. .

{\begin{matrix} ϕ^{k} (x, t) = 0, x \in Γ^{k} (t), k = V, M, folding lines \\ ϕ^{k} (x, t) \neq 0, x \in D ∖ Γ^{k} (t), facet \end{matrix}

(1)

where 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. , which is defined by differentiating the level set function with respect to time t [44].

\frac{\partial ϕ^{k}}{\partial t} - V^{k} \cdot \nabla ϕ^{k} = 0

(2)

where

V^{k} = \dot{x}

is the velocity field for the folding line. During the topology optimization, each level set function involves individually. Thus, in the rest of the article, for simplicity, we use φ to represent a general level set function unless otherwise stated.

Fig. 2

A schematic of the level set representation of an origami structure

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.

The origami mechanism is simulated using shell models where the in-plane stretching, out-of-plane bending, and shear deformation can be well captured. For simplicity, at this point, we assume the problem to be linear elastic. The crease is considered a region of locally weak material with directional imperfections on the shell. A similar technique has been used in studying shells with curved creases [28] and the collapse of a thin-walled tubular [49]. We define the thickness and Young’s modules on the crease as half of those in the facet area.

\begin{aligned} t_{f a c t e s} = 2 t_{f o l d s} \\ E_{f a c t e s} = 2 E_{f o l d s} \end{aligned}

(3)

where 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].

δ (ϕ) = {\begin{matrix} \frac{1}{2 ε} + \frac{1}{2 ε} \cos (\frac{π ϕ}{ε}) - ε \leq ϕ \leq ε \\ 0 else \end{matrix}

(4)

where ɛ controls the size of the transition area, which is chosen as two times the mesh size [47,43]. Then, we normalize δ(φ) to

\hat{δ} (ϕ)

⁠, where the maximum value is equal to 1, and couple the

\hat{δ} (ϕ)

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 (ξ^{1}, ξ^{2}) = r_{R} (ξ^{1}, ξ^{2}) + ζ_{0} n (ξ^{1}, ξ^{2})

(5)

where r_R is the meshed surface, ξ₁ and ξ₂ are the two curvilinear coordinates in the middle surface, and n is the normal direction on the shell. The ζ₀ is the offset in the thickness direction, which is related to

\hat{δ} (ϕ)

as shown in Eq. .

ζ_{0} = {\begin{matrix} + \hat{δ} (ϕ)_{M} \times \frac{1}{2} \times t_{f o l d s}, on mountain folds \\ 0, on facet area \\ - \hat{δ} (ϕ)_{V} \times \frac{1}{2} \times t_{f o l d s}, on valley folds \end{matrix}

(6)

Fig. 3

The schematic of level set representation of a modified shell model

The strain of the undeformed middle surface can be evaluated as follows:

\frac{\partial r}{\partial ξ^{α}} = \frac{\partial r_{R}}{\partial ξ^{α}} + ζ_{0} \frac{\partial n}{\partial {\bar{ξ}}^{α}}

(7)

where α 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 (ζ₀) have been canceled out. Thus, the overlaps are essentially areas with weak material but without the in-thickness offsets.

Fig. 4

Origami structures modeled by modified shell: (a) modeling of origami structures with straight folds and (b) modelling of origami structures with curved folds

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₁) to meet the kinematic requirements;
the characteristic requirements for the design to be an origami structure (J₂);
the perimeter constraint for achieving concise designs.

The general objective function can be formulated as follows:

\begin{aligned} Minimize : J = I_{1} J_{1} + I_{2} J_{2} + α | \partial Ω |, \\ Subject to : a (u, v) = l (v), \forall v \in U_{a d} \end{aligned}

(8)

In the aforementioned equation, U_ad is the space kinematically admissible displacement [55], Ω is the region occupied with the linear elastic material. I₁, I₂, and α are constant weighting factors for kinematic target, characteristic objective, and perimeter constraint, respectively. The internal virtual energy a(u, v) and the exterior virtual energy l(v) are defined as follows:

\begin{aligned} a (u, v) = \int_{Ω} ε_{i j} (u) C_{i j k l} ε_{k l} (v) d Ω, \\ l (v) = \int_{Ω} f \cdot v d Ω + \int_{Γ_{N}} g \cdot v d Γ \end{aligned}

(9)

where D is the design domain; Γ = ∂Ω is the boundary of the design; a(u, v) is a symmetric bilinear function in terms of the displacement u, and the test function v, that is, a(u, v) = a(v, u); and l is a linear function in terms of the body force f and the traction force g on the Neumann boundary conditions Γ_N, 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 κ and can be calculated by [56]:

κ = - \nabla \cdot (\frac{\nabla ϕ}{| \nabla ϕ |})

(10)

Fig. 5

A schematic of general boundary condition

Therefore, we will only show the detailed steps for sensitivity analysis of the kinematic (J₁) and the characteristic requirements (J₂).

3.1 Shape Sensitivity Analysis of the Kinematic Objective.

Here, we adopt the classic objective function of kinematic requirement–the least square minimization between the displacement u and a given target displacement u₀ of a selected zone w(x) in the design domain; w(x) is zero excepted on the selected area.

\begin{aligned} Minimize : J_{1} = {[\int_{Ω} w (x) {| u - u_{0} |}^{2} d Ω]}^{\frac{1}{2}} \\ Subject to : a (u, v) = l (v), \forall v \in U_{a d} \end{aligned}

(11)

The Lagrangian of the optimization problem is expressed as follows:

L_{1} (u, v) = J_{1} + a (u, v) - l (v)

(12)

The material derivative of the Lagrangian equation with respect to a pseudo time t is formulated as follows:

\frac{d L_{1} (u, v)}{d t} = \frac{\partial L_{1} (u, v)}{\partial t} + \frac{\partial L_{1} (u, v)}{\partial Ω}

(13)

where the partial derivative with respect to time results in the so-called adjoint equation:

\frac{\partial L_{1} (v)}{\partial t} = L_{1}^{^{'}} = c_{0} \int_{Ω} w (x) (u - u_{0}) \cdot u^{'} d Ω + a (u^{'}, v)

(14)

where

c_{0} = {[\int_{Ω} w (x) {| u - u_{0} |}^{2} d Ω]}^{- \frac{1}{2}}

(15)

The convection term of the material derivative forms the shape derivative:

\begin{aligned} \frac{\partial L_{1} (u, v)}{\partial Ω} & = \int_{Γ} \frac{c_{0}}{2} w (x) {| u - u_{0} |}^{2} V_{n 1} d s \\ + \int_{Γ} ε_{i j} (u) C_{i j k l} ε_{k l} (v) V_{n 1} d s - \int_{Γ} f \cdot v V_{n 1} d s \\ - \int_{Γ} [\frac{\partial (g \cdot v)}{\partial n} + κ g \cdot v] V_{n 1} d s \end{aligned}

(16)

Using the steepest descent method while neglecting body force and traction force, the normal velocity can be written as follows:

V_{n 1} = - [\frac{c_{0}}{2} w (x) {| u - u_{0} |}^{2} + ε_{i j} (u) C_{i j k l} ε_{k l} (v)]

(17)

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.

The J_2I and J_2II can be represented as following equations:

\begin{matrix} J_{2_{I}} = a (u, u) - a_{Γ} (u, u), \\ J_{2_{II}} = \frac{1}{a_{b} (u, u)} \end{matrix}

(18)

where a(u, u) is the total strain energy, a_Γ(u, u) is the total strain energy on the creases, and a_b(u, u) is the total bending energy.

\begin{aligned} a_{Γ} (u, v) = \int_{τ} ε_{i j} (u) C_{i j k l} ε_{k l} (u) d s, \\ a_{b} (u, u) = \int_{Ω} ε_{i j}^{b} (u) C_{i j k l} ε_{k l}^{b} (u) d Ω \end{aligned}

(19)

We can formulate the characteristic objective function J₂ as follows:

J_{2} = J_{2_{I}} J_{2_{II}}

(20)

Then, the shape derivative of the J₂ can be represented as follows:

D_{Ω} J_{2} = J_{2_{II}} D_{Ω} J_{2_{I}} + J_{2_{I}} D_{Ω} J_{2_{II}}

(21)

where 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.

We summarize the first condition as an objective function J_2I, which minimizes the difference of total strain energy a(u, u) and the total strain energy on the creases a_Γ(u, u), as shown in Eq. :

\begin{aligned} Minimize : = J_{2_{I}} a (u, u) - a_{Γ} (u, u), \\ Subject to : a (u, v) = l (v), \forall v \in U_{a d} \end{aligned}

(22)

The Lagrangian of the optimization problem can be formulated as follows:

L_{2} = a (u, v) - a_{Γ} (u, v) + a (u, v) - l (v)

(23)

The material time derivative of the Lagrangian equation can be calculated as follows:

\frac{d L_{2}}{d t} = \frac{\partial L_{2}}{\partial t} + \frac{\partial L_{2}}{\partial Ω}

(24)

\frac{\partial L_{2}}{\partial t} = 2 a (u^{'}, u) - 2 a_{Γ} (u^{'}, u) + a (u^{'}, v)

(25)

\begin{aligned} \frac{\partial L_{2}}{\partial Ω} & = \int_{τ} ε_{i j} (u) C_{i j k l} ε_{k l} (u) V_{n 2 I} d s - \int_{τ} [\frac{\partial ε_{i j} (u) C_{i j k l} ε_{k l} (u)}{\partial n} + κ a (u, u)] \\ \times V_{n 2 I} d s + \int_{τ} ε_{i j} (u) C_{i j k l} ε_{k l} (v) V_{n 2 I} d s \end{aligned}

(26)

The adjoint equation Eq. can be expressed as follows:

\begin{aligned} \frac{\partial L_{2}}{\partial t} & = 2 \int_{Ω} ε_{i j} (u^{'}) C_{i j k l} ε_{k l} (u) d Ω - 2 \int_{Ω} ε_{i j} (u^{'}) C_{i j k l} ε_{k l} (u) δ (ϕ) \\ \times | \nabla ϕ (x) | d Ω + \int_{Ω} ε_{i j} (u^{'}) C_{i j k l} ε_{k l} (v) d Ω \end{aligned}

(27)

where v is the adjoint variable and is achieved by solving Eq. using the finite element method (FEM).

3.2.2 Maximizing the Out-of-Plane Deformation.

The second characteristic requirement 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:

\begin{aligned} Minimize : J_{2_{II}} = \frac{1}{a_{b} (u, u)}, \\ Subject to : a (u, v) = l (v), \forall v \in U_{a d} \end{aligned}

(28)

The optimization problem becomes:

L_{3} = \frac{1}{a_{b} (u, u)} + a (u, v) - l (v)

(29)

The material time derivative can be formulated as follows:

\frac{d L_{3}}{d t} = \frac{\partial L_{3}}{\partial t} + \frac{\partial L_{3}}{\partial Ω}

(30)

where

\begin{matrix} \frac{\partial L_{3}}{\partial t} = \frac{- 2 a_{b} (u^{'}, u)}{a_{b} {(u, u)}^{2}} + a (u^{'}, v) - β \int_{Γ} u^{'} v d Γ_{D}, \\ \frac{\partial L_{3}}{\partial Ω} = \frac{- 1}{a^{2} (u, u)} \int_{Γ} ε_{i j}^{b} (u) C_{i j k l} ε_{k l}^{b} (u) V_{n 2 I I} d s + \int_{Γ} ε_{i j} (u) C_{i j k l} ε_{k l} (v) V_{n 2 I I} d s \end{matrix}

(31)

To solve the adjoint variable v, we solve the adjoint equation:

a (u^{'}, v) - \frac{2 a_{b} (u^{'}, u)}{a_{b} {(u, u)}^{2}} - β \int_{Γ} u^{'} v d Γ_{D} = 0

(32)

Now we can substitute the shape gradient of each section into Eq. , assuming body force is zero and ignoring the velocity along the zero-length Neumann boundary, the shape gradient of the objective function is expressed as follows:

D_{Ω} J_{2} = \int_{Γ} G V_{n 2} d Γ_{D}

(33)

Thus, by applying the steepest descent method, the normal velocity can be written as follows:

V_{n 2} = - G

(34)

where

\begin{aligned} G = a_{b} (u, u) [ε_{i j} (u) C_{i j k l} ε_{k l} (u) - \frac{\partial ε_{i j} (u) C_{i j k l} ε_{k l} (u)}{\partial n} - κ a (u, u) + ε_{i j} (u) C_{i j k l} ε_{k l} (v)] \\ + [a (u, u) - a_{Γ} (u, u)] [\frac{- 1}{a^{2} (u, u)} ε_{i j}^{b} (u) C_{i j k l} ε_{k l}^{b} (u) + ε_{i j} (u) C_{i j k l} ε_{k l} (v)] \end{aligned}

(35)

4 Numerical Implementation

4.1 An Origami Gripper Design.

In this section, we present an optimization example to test the method's performance. The objective is to obtain an origami mechanism with a minimum crease number where the target area(s) would reach a displacement in the target direction. In other words, the optimized structure can transfer force and displacement from the input port area to the output port by deformation on the optimized origami.

\begin{aligned} Minimize : J = I_{1} {[\int_{Ω} w (x) {| u - u_{0} |}^{2} d Ω]}^{\frac{1}{2}} \\ + I_{2} \frac{a (u, u) - a_{Γ} (u, u)}{a_{b} (u, u)} + α | \partial Ω |, \\ Subject to : a (u, v) = l (v), \forall v \in U_{a d} \end{aligned}

(36)

where I₁, I₂, 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₂ = 1−I₂. Thus, the velocity field can be represented as follows:

V_{n} = I_{1} V_{n 1} + I_{2} V_{n 2} + α κ

(37)

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₀ 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⁶ Pa on the facet and 5 × 10⁵ Pa on the folding lines.

Fig. 6

Topology optimization of an origami twister: (a) boundary conditions, (b) the evolution on 2D, and (c) the deformed status on 3D

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.

Fig. 7

The optimization history of the origami gripper

Figure 8(a) shows the fabricated prototype. As shown in Fig. 8(b), the top and bottom of the gripper can be closed up by pushing the middle portion of both sides.

Fig. 8

Fabricated prototype: (a) the fabricated origami gripper and (b) the folded origami gripper

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⁶ Pa on the facet and 5 × 10⁵ Pa on the folding lines.

Fig. 9

Topology optimization of an origami twister: (a) boundary conditions, (b) the evolution on 2D, and (c) the deformed status on 3D

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.

Fig. 10

Prototype of the twister: (a) the fabricated twister and (b) the folded status

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.

Fig. 11

Design evolution with different initializations: (a) less initial folds and (b) packed initial folds

5.2 The Performance of the Characteristic Objective.

As we discussed in Section 3.2, the characteristic objective is an energy function essential for reaching a design where most deformation happens on the folding lines. An example is presented to show the effect of the characteristic objective J₂. The objective function only contains the J₂ with the linear elastic constraint, as shown in the equations:

\begin{aligned} Minimize : J_{2} = \int_{Ω} \frac{a (u, u) - a_{Γ} (u, u)}{a_{b} (u, u)} \\ Subject to : a (u, v) = l (v), \forall v \in U_{ad} \end{aligned}

(38)

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.

Fig. 12

(a) Boundary conditions, (b) the initial design, and (c) the final design

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.

Fig. 13

Left to right: (a) convergence plots and (b) strain energy density field on the final design

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

In this section, the background information of the linear elastic shell structure is introduced, and the derivation of the strain energy density is presented. Typically, a shell structure with a curved external surface and a uniform thickness is shown in Fig. 14. Let ξ₁ and ξ₂ be the two curvilinear coordinates in the middle surface and ξ₃ 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(ξ₁, ξ₂, ξ₃) [63,64].

x (ξ_{1}, ξ_{2}, ξ_{3}) = x^{n} (ξ_{1}, ξ_{2}) + ξ_{3} \frac{t (ξ_{1}, ξ_{2})}{2} n (ξ_{1}, ξ_{2})

(A1)

Fig. 14

The schematic figure of shell structure

The Jacobian matrix of the mapping between the middle surface and reference surface is shown in the left corner in Fig. 14, where

J = \frac{\partial x_{i}}{\partial ξ_{j}}

(A2)

For linear elastic material, the strain tensor is defined as follows:

\begin{matrix} ε_{i j} (u) = \frac{1}{2} (\frac{\partial u_{i}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial x_{i}}) \\ ε_{i j} (u) = \frac{1}{2} (\frac{\partial u_{i}}{\partial ξ_{m}} \frac{\partial ξ_{m}}{\partial x_{j}} + \frac{\partial u_{j}}{\partial ξ_{m}} \frac{\partial ξ_{m}}{\partial x_{i}}) \end{matrix}

(A3)

\begin{matrix} ε_{i j} (u) = \frac{1}{2} (\frac{\partial u_{i}}{\partial ξ_{m}} \frac{1}{J_{m j}} + \frac{\partial u_{j}}{\partial ξ_{m}} \frac{1}{J_{m i}}) \\ = sym (\frac{\partial u_{i}}{\partial ξ_{m}} J_{m j}^{- 1}) \end{matrix}

(A4)

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

a (u, v) = \int_{Ω} ε_{i j} (u) C_{i j k l} ε_{k l} (v) d Ω

(A4)

where u stands for the displacement and v is the test function.

According to the Reissner–Mindlin theory of plate [65], we can assume that the displacement along the thickness is linear, and the cross section of the shell structure remains flat after deformation. Then, the deformation on the shell can be written as follows [66]:

u = u^{1} (ξ_{1}, ξ_{2}) + ξ_{3} u^{2} (ξ_{1}, ξ_{2})

(A5)

where on the right-hand side, the first term represents the membrane deformation, and the second term stands for the bending and shear deformation.

For the thin shell, whose thickness is very small and can be neglected, the Jacobian is considered to be a function of only ξ₁,ξ₂ coordinates. Thus, Eq. (42) can be simplified as follows [66]:

ε_{i j} (u) = ε_{i j}^{1} (u) + ξ_{3} ε_{i j}^{2} (u)

(A6)

where ɛ_ij¹(u) and ɛ_ij²(u) are the membrane-shear stain and the bending strain, respectively.

Substitute Eq. (45) into (43) and integrate from ξ₃ ∈ [−1, 1], we can rewrite the energy bilinear form as follows [55,66]:

\begin{aligned} a (u, v) & = \int_{- 1}^{1} \int \int [ε_{i j}^{1} (u) + ξ_{3} ε_{i j}^{2} (u)] C_{i j k l} [ε_{i j}^{1} (v) + ξ_{3} ε_{i j}^{2} (v)] | J | d ξ_{1} d ξ_{2} d ξ_{3} \\ = \int_{- 1}^{1} \int \int [ε_{i j}^{1} (u) C_{i j k l} ε_{i j}^{1} (v) + ξ_{3} ε_{i j}^{1} (u) C_{i j k l} ε_{i j}^{2} (v) \\ + ξ_{3} ε_{i j}^{2} (u) C_{i j k l} ε_{i j}^{1} (v) + ξ_{3}^{2} ε_{i j}^{2} (u) C_{i j k l} ε_{i j}^{2} (v)] | J | d ξ_{1} d ξ_{2} d ξ_{3} \end{aligned}

(A7)

since the second and third terms are odd functions over the interval ξ₃ ∈ [−1, 1], and thus, Eq. (46) can be simplified as follows [55]:

a (u, v) = \int_{A} [2 ε_{i j}^{1} (u) C_{i j k l} ε_{k l}^{1} (v) | J | + \frac{2}{3} ε_{i j}^{2} (u) C_{i j k l} ε_{k l}^{2} (v) | J |] d A

(A8)

where the first term represents the coupled membrane-shear strain energy and the second term donates the bending strain energy. We use a_b(u, u) to denote the bending strain energy in the contents mentioned earlier.

References

1.

Turner

,

N.

,

Goodwine

,

B.

, and

Sen

,

M.

,

2016

, “

A Review of Origami Applications in Mechanical Engineering

,”

Proc. Inst. Mech. Eng., Part C

,

230

(

14

), pp.

2345

–

2362

.

2.

Overvelde

,

J. T.

,

De Jong

,

T. A.

,

Shevchenko

,

Y.

,

Becerra

,

S. A.

,

Whitesides

,

G. M.

,

Weaver

,

J. C.

,

Hoberman

,

C.

, and

Bertoldi

,

K.

,

2016

, “

A Three-Dimensional Actuated Origami-Inspired Transformable Metamaterial With Multiple Degrees of Freedom

,”

Nat. Commun.

,

7

(

1

), pp.

1

–

8

.

3.

Fuchi

,

K.

,

Diaz

,

A. R.

,

Rothwell

,

E. J.

,

Ouedraogo

,

R. O.

, and

Tang

,

J.

,

2012

, “

An Origami Tunable Metamaterial

,”

J. Appl. Phys.

,

111

(

8

), p.

084905

.

4.

Kamrava

,

S.

,

Mousanezhad

,

D.

,

Ebrahimi

,

H.

,

Ghosh

,

R.

, and

Vaziri

,

A.

,

2017

, “

Origami-Based Cellular Metamaterial With Auxetic, Bistable, and Self-Locking Properties

,”

Sci. Rep.

,

7

(

1

), pp.

1

–

9

.

PubMed

5.

Morgan

,

J.

,

Magleby

,

S. P.

, and

Howell

,

L. L.

,

2016

, “

An Approach to Designing Origami-Adapted Aerospace Mechanisms

,”

ASME J. Mech. Des.

,

138

(

5

), p.

052301

.

6.

Ishida

,

S.

,

Nojima

,

T.

, and

Hagiwara

,

I.

,

2014

, “

Mathematical Approach to Model Foldable Conical Structures Using Conformal Mapping

,”

ASME J. Mech. Des.

,

136

(

9

), p.

091007

.

7.

Hull

,

T.

,

1994

, “

On the Mathematics of Flat Origamis

,”

Congressus Numerantium

,

Boca Raton, FL

,

Mar. 7–11

, pp.

215

–

224

.

8.

Hull

,

T. C.

,

2002

, “

Modelling the Folding of Paper Into Three Dimensions Using Affine Transformations

,”

Linear Algebra Appl.

,

348

(

1–3

), pp.

273

–

282

.

9.

Lang

,

R. J.

,

1996

, “

A Computational Algorithm for Origami Design

,”

Proceedings of the Twelfth Annual Symposium on Computational Geometry

,

Philadelphia, PA

,

May 24–26

, ACM, pp.

98

–

105

.

10.

Tachi

,

T.

,

2010

, “

Geometric Considerations for the Design of Rigid Origami Structures

,”

Proceedings of the International Association for Shell and Spatial Structures (IASS) Symposium

,

Shanghai, China

,

Nov. 8–12

, Vol. 12, pp.

458

–

460

.

11.

Schenk

,

M.

, and

Guest

,

S. D.

,

2010

, “

Origami Folding: A Structural Engineering Approach

,”

Origami 5: Fifth International Meeting of Origami Science, Mathematics, and Education

,

Singapore

,

July 13–17

, CRC Press, Boca Raton, FL, pp.

291

–

304

.

12.

Tachi

,

T.

,

2009

, “Simulation of Rigid Origami,”

Origami 4

, Vol.

4

,

CRC Press

,

New York

, pp.

175

–

187

.

13.

Wei

,

Z. Y.

,

Guo

,

Z. V.

,

Dudte

,

L.

,

Liang

,

H. Y.

, and

Mahadevan

,

L.

,

2013

, “

Geometric Mechanics of Periodic Pleated Origami

,”

Phys. Rev. Lett.

,

110

(

21

), p.

215501

.

PubMed

14.

Lang

,

R. J.

,

Tolman

,

K. A.

,

Crampton

,

E. B.

,

Magleby

,

S. P.

, and

Howell

,

L. L.

,

2018

, “

A Review of Thicknessaccommodation Techniques in Origami-Inspired Engineering

,”

ASME Appl. Mech. Rev.

,

70

(

1

), p.

010805

.

15.

Peraza Hernandez

,

E. A.

,

Hartl

,

D. J.

, and

Lagoudas

,

D. C.

,

2016

, “

Kinematics of Origami Structures with Smooth Folds

,”

ASME J. Mech. Rob.

,

8

(

6

), p.

061019

.

16.

Hernandez

,

E. A. P.

,

Hartl

,

D. J.

, and

Lagoudas

,

D. C.

,

2019

,

Active Origami

,

Springer Cham

,

Switzerland

, pp.

201

–

268

.

17.

Liu

,

K.

, and

Paulino

,

G. H.

,

2016

, “

Merlin: A Matlab Implementation to Capture Highly Nonlinear Behavior of Nonrigid Origami

,”

Proceedings of IASS Annual Symposia

,

Tokyo, Japan

,

Sept. 26

, pp.

1

–

10

.

18.

Liu

,

K.

, and

Paulino

,

G.

,

2017

, “

Nonlinear Mechanics of non-Rigid Origami: An Efficient Computational Approach

,”

Proc. R. Soc. A

,

473

(

2206

), p.

20170348

.

19.

Filipov

,

E.

,

Liu

,

K.

,

Tachi

,

T.

,

Schenk

,

M.

, and

Paulino

,

G.

,

2017

, “

Bar and Hinge Models for Scalable Analysis of Origami

,”

Int. J. Solids Struct.

,

124

, pp.

26

–

45

.

20.

Ghassaei

,

A.

,

Demaine

,

E. D.

, and

Gershenfeld

,

N.

,

2018

, “

Fast, Interactive Origami Simulation Using gpu Computation

,”

7th International Meeting on Origami in Science, Mathematics and Education

,

Oxford, UK

,

Sept. 4–7

.

21.

Gilewski

,

W.

,

Pełczynski

,

J.

, and

Stawarz

,

P.

,

2014

, “

A Comparative Study of Origami Inspired Folded Plates

,”

Procedia Eng.

,

91

, pp.

220

–

225

.

22.

Cai

,

J.

,

Ren

,

Z.

,

Ding

,

Y.

,

Deng

,

X.

,

Xu

,

Y.

, and

Feng

,

J.

,

2017

, “

Deployment Simulation of Foldable Origami Membrane Structures

,”

Aerosp. Sci. Technol.

,

67

, pp.

343

–

353

.

23.

Dias

,

M. A.

,

Dudte

,

L. H.

,

Mahadevan

,

L.

, and

Santangelo

,

C. D.

,

2012

, “

Geometric Mechanics of Curved Crease Origami

,”

Phys. Rev. Lett.

,

109

(

11

), p.

114301

.

PubMed

24.

Kilian

,

M.

,

Flory

,

S.

,

Chen

,

Z.

,

Mitra

,

N. J.

,

Sheffer

,

A.

, and

Pottmann

,

H.

,

2008

, “

Curved Folding

,”

ACM Trans. Graph.

,

27

(

3

), pp.

1

–

9

.

25.

Kilian

,

M.

,

Monszpart

,

A.

, and

Mitra

,

N. J.

,

2017

, “

String Actuated Curved Folded Surfaces

,”

ACM Trans. Graph.

,

36

(

3

), pp.

1

–

13

.

26.

Vergauwen

,

A.

,

De Laet

,

L.

, and

De Temmerman

,

N.

,

2017

, “

Computational Modelling Methods for Pliable Structures Based on Curved-Line Folding

,”

Comput. Aided Des.

,

83

(

C

), pp.

51

–

63

.

27.

Callens

,

S. J.

, and

Zadpoor

,

A. A.

,

2018

, “

From Flat Sheets to Curved Geometries: Origami and Kirigami Approaches

,”

Mater. Today

,

21

(

3

), pp.

241

–

264

.

28.

Bende

,

N. P.

,

Evans

,

A. A.

,

Innes-Gold

,

S.

,

Marin

,

L. A.

,

Cohen

,

I.

,

Hayward

,

R. C.

, and

Santangelo

,

C. D.

,

2015

, “

Geometrically Controlled Snapping Transitions in Shells With Curved Creases

,”

Proc. Natl. Acad. Sci. USA

,

112

(

36

), pp.

11175

–

11180

.

29.

Nelson

,

T. G.

,

Lang

,

R. J.

,

Magleby

,

S. P.

, and

Howell

,

L. L.

,

2016

, “

Curved-Folding-Inspired Deployable Compliant Rolling-Contact Element (d-Core)

,”

Mech. Mach. Theory

,

96

(

Part 2

), pp.

225

–

238

.

30.

Badger

,

J. C.

,

Nelson

,

T. G.

,

Lang

,

R. J.

,

Halverson

,

D. M.

, and

Howell

,

L. L.

,

2019

, “

Normalized Coordinate Equations and an Energy Method for Predicting Natural Curved-Fold Configurations

,”

ASME J. Appl. Mech.

,

86

(

7

), p.

071006

.

31.

Greenwood

,

J.

,

Avila

,

A.

,

Howell

,

L.

, and

Magleby

,

S.

,

2020

, “

Conceptualizing Stable States in Origami-Based Devices Using an Energy Visualization Approach

,”

ASME J. Mech. Des.

,

142

(

9

), p.

093302

.

32.

Avila

,

A.

,

Greenwood

,

J.

,

Magleby

,

S. P.

, and

Howell

,

L. L.

,

2019

, “

Conceptualizing Stable States in Origami-Based Devices Using an Energy Visualization Approach

,”

ASME 2019 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference

,

Anaheim, CA

,

Aug. 18–21

.

33.

Fuchi

,

K.

, and

Diaz

,

A. R.

,

2013

, “

Origami Design by Topology Optimization

,”

ASME J. Mech. Des.

,

135

(

11

), p.

111003

.

34.

Fuchi

,

K.

,

Buskohl

,

P. R.

,

Bazzan

,

G.

,

Durstock

,

M. F.

,

Reich

,

G. W.

,

Vaia

,

R. A.

, and

Joo

,

J. J.

,

2015

, “

Origami Actuator Design and Networking Through Crease Topology Optimization

,”

ASME J. Mech. Des.

,

137

(

9

), p.

091401

.

35.

Fuchi

,

K.

,

Buskohl

,

P. R.

,

Bazzan

,

G.

,

Durstock

,

M. F.

,

Reich

,

G. W.

,

Vaia

,

R. A.

, and

Joo

,

J. J.

,

2016

, “

Design Optimization Challenges of Origami-Based Mechanisms With Sequenced Folding

,”

ASME J. Mech. Rob.

,

8

(

5

), p.

051011

.

36.

Gillman

,

A. S.

,

Fuchi

,

K.

, and

Buskohl

,

P. R.

,

2019

, “

Discovering Sequenced Origami Folding Through Nonlinear Mechanics and Topology Optimization

,”

ASME J. Mech. Des.

,

141

(

4

), p.

041401

.

37.

Yu

,

Y.

,

Hong

,

T.-C. K.

,

Economou

,

A.

, and

Paulino

,

G. H.

,

2021

, “

Rethinking Origami: A Generative Specification of Origami Patterns With Shape Grammars

,”

Comput. Aided Des.

,

137

, p.

103029

.

38.

Suzuki

,

K.

, and

Kikuchi

,

N.

,

1991

, “

A Homogenization Method for Shape and Topology Optimization

,”

Comput. Methods Appl. Mech. Eng.

,

93

(

3

), pp.

291

–

318

.

39.

Allaire

,

G.

, and

Kohn

,

R.

,

1993

, “Topology Optimization and Optimal Shape Design Using Homogenization,”

Topology Design of Structures

,

M.P.

Bendsøe

, and

C. A. M.

Soares

, eds.,

Springer

,

Dordrecht

, pp.

207

–

218

.

40.

Nishiwaki

,

S.

,

Frecker

,

M. I.

,

Min

,

S.

, and

Kikuchi

,

N.

,

1998

, “

Topology Optimization of Compliant Mechanisms Using the Homogenization Method

,”

Int. J. Numer. Methods Eng.

,

42

(

3

), pp.

535

–

559

.

41.

Rozvany

,

G. I.

,

Zhou

,

M.

, and

Birker

,

T.

,

1992

, “

Generalized Shape Optimization Without Homogenization

,”

Struct. Optim

,

4

(

3–4

), pp.

250

–

252

.

42.

Bendsoe

,

M. P.

, and

Sigmund

,

O.

,

2013

,

Topology Optimization: Theory, Methods, and Applications

,

Springer Science & Business Media

,

Germany

.

43.

van Dijk

,

N. P.

,

Maute

,

K.

,

Langelaar

,

M.

, and

Van Keulen

,

F.

,

2013

, “

Level-Set Methods for Structural Topology Optimization: A Review

,”

Struct. Multidiscipl. Optim.

,

48

(

3

), pp.

437

–

472

.

44.

Allaire

,

G.

,

Jouve

,

F.

, and

Toader

,

A.-M.

,

2004

, “

Structural Optimization Using Sensitivity Analysis and a Level-Set Method

,”

J. Comput. Phys.

,

194

(

1

), pp.

363

–

393

.

45.

Osher

,

S. J.

, and

Santosa

,

F.

,

2001

, “

Level Set Methods for Optimization Problems Involving Geometry and Constraints: I. Frequencies of a Two-Density Inhomogeneous Drum

,”

J. Comput. Phys.

,

171

(

1

), pp.

272

–

288

.

46.

Allaire

,

G.

,

Jouve

,

F.

, and

Toader

,

A.-M.

,

2002

, “

A Level-Set Method for Shape Optimization

,”

Comptes Rendus Math.

,

334

(

12

), pp.

1125

–

1130

.

47.

Wang

,

M. Y.

,

Wang

,

X.

, and

Guo

,

D.

,

2003

, “

A Level Set Method for Structural Topology Optimization

,”

Comput. Methods Appl. Mech. Eng.

,

192

(

1–2

), pp.

227

–

246

.

48.

Wang

,

M. Y.

, and

Wang

,

S.

,

2006

, “Parametric Shape and Topology Optimization With Radial Basis Functions,”

IUTAM Symposium on Topological Design Optimization of Structures, Machines and Materials

,

Springer

, pp.

13

–

22

.

49.

Bandi

,

P.

,

Detwiler

,

D.

,

Schmiedeler

,

J. P.

, and

Tovar

,

A.

,

2015

, “

Design of Progressively Folding Thin-Walled Tubular Components Using Compliant Mechanism Synthesis

,”

Thin-Walled Struct.

,

95

, pp.

208

–

220

.

50.

Zhu

,

J.

,

2013

,

The Finite Element Method: Its Basis and Fundamentals

,

Elsevier

,

New York

.

51.

Zienkiewicz

,

O. C.

, and

Taylor

,

R. L.

,

2005

,

The Finite Element Method for Solid and Structural Mechanics

,

Elsevier

,

New York

.

52.

Liu

,

G.-R.

, and

Quek

,

S. S.

,

2013

,

The Finite Element Method: A Practical Course

,

Elsevier Science

,

Netherlands

.

53.

Bathe

,

K.-J.

, and

Chapelle

,

D.

,

2003

,

Finite Element Analysis of Shells–Fundamentals

,

Springer-Verlag Telos

,

Germany

.

54.

Reddy

,

J. N.

,

2003

,

Mechanics of Laminated Composite Plates and Shells: Theory and Analysis

,

CRC Press

,

Boca Raton, FL

.

55.

Choi

,

K. K.

, and

Kim

,

N.-H.

,

2006

,

Structural Sensitivity Analysis and Optimization 1: Linear Systems

,

Springer Science & Business Media

,

New York

.

56.

Osher

,

S.

, and

Fedkiw

,

R.

,

2006

,

Level set Methods and Dynamic Implicit Surfaces

, Vol.

153

,

Springer Science & Business Media

,

New York

.

57.

Tian

,

J.

,

Li

,

M.

,

Han

,

Z.

,

Chen

,

Y.

,

Gu

,

X. D.

,

Ge

,

Q.

, and

Chen

,

S.

,

2022

, “

Conformal Topology Optimization of Multi-Material Ferromagnetic Soft Active Structures Using an Extended Level Set Method

,”

Comput. Methods Appl. Mech. Eng.

,

389

, p.

114394

.

58.

Ye

,

Q.

,

Guo

,

Y.

,

Chen

,

S.

,

Gu

,

X. D.

, and

Lei

,

N.

,

2018

, “

Topology Optimization of Conformal Structures Using Extended Level Set Methods and Conformal Geometry Theory

,”

ASME 2018 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference

,

Quebec City, Canada

,

Aug. 26–29

.

59.

Ye

,

Q.

,

Guo

,

Y.

,

Chen

,

S.

,

Lei

,

N.

, and

Gu

,

X. D.

,

2019

, “

Topology Optimization of Conformal Structures on Manifolds Using Extended Level Set Methods (x-lsm) and Conformal Geometry Theory

,”

Comput. Methods Appl. Mech. Eng.

,

344

, pp.

164

–

185

.

60.

Ye

,

Q.

,

Jiang

,

L.

,

Chen

,

S.

, and

Gu

,

X.

,

2018

, “

Generative Design of Multifunctional Conformal Structures Using Extended Level Set Methods (x-lsm) and Conformal Geometry Theory

,”

IUTAM Symposium on When Topology Optimization Meets Additive Manufacturing—Theory and Methods.

,

Dalian, China

,

Oct. 8–12

.

61.

Ye

,

Q.

,

Gu

,

X. D.

, and

Chen

,

S.

,

2020

, “

Generative Design of Origami-Inspired Mechanisms With a Variational Level Set Approach

,”

ASME 2020 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference

,

Online, Virtual

,

Aug. 17–19

.

62.

Xu

,

X.

,

Chen

,

S.

,

Gu

,

X. D.

, and

Wang

,

M. Y.

,

2021

, “

Conformal Topology Optimization of Heat Conduction Problems on Manifolds Using an Extended Level Set Method (x-lsm)

,”

International Design Engineering Technical Conferences and Computers and Information in Engineering Conference

,

Virtual, Online

,

Aug. 17–19

.

63.

Ciarlet

,

P. G.

,

2005

, “

An Introduction to Differential Geometry With Applications to Elasticity

,”

J. Elast.

,

78

(

1–3

), pp.

1

–

215

.

64.

Chapelle

,

D.

, and

Bathe

,

K.-J.

,

2010

,

The Finite Element Analysis of Shells-Fundamentals

,

Springer Science & Business Media

,

Berlin/Heidelberg

.

65.

Papadopoulos

,

P.

, and

Taylor

,

R. L.

,

1990

, “

A Triangular Element Based on Reissner-Mindlin Plate Theory

,”

Int. J. Numer. Methods Eng.

,

30

(

5

), pp.

1029

–

1049

.

66.

Kim

,

N. H.

,

Choi

,

K. K.

,

Chen

,

J.-S.

, and

Botkin

,

M. E.

,

2002

, “

Meshfree Analysis and Design Sensitivity Analysis for Shell Structures

,”

Int. J. Numer. Methods Eng.

,

53

(

9

), pp.

2087

–

2116

.