## Abstract

Active origami structures usually have creases made with soft and compliant plates because it is difficult to fabricate real hinges and actuate them. However, most conventional origami modeling techniques fail to capture these compliant creases and simplify them as concentrated rotational springs, which neglects torsional and extensional deformations of the creases. In this paper, an improved formulation of a bar and hinge model is proposed to explicitly capture the geometry and the flexibility of compliant creases with nonnegligible width in an origami, and the model is verified against finite element simulations. The verification shows that the model performs relatively well despite being simple and computationally inexpensive. Moreover, simulation examples demonstrate that the proposed model can capture the bistable behavior of the compliant crease origami with nonnegligible crease width because it explicitly includes the extensional stretching energy into the simulation framework and allows torsional crease deformations.

## Introduction

Origami-inspired assemblages provide feasible and convenient approaches for designing reconfigurable structures and metamaterials with tunable mechanical properties. Various engineering structures of this type have been proposed in the past decades, such as reconfigurable facades and canopies [1], bistable wings [2], mechanical memory devices [3], pneumatic actuators [4], and various metamaterials [5,6]. To understand the mechanical and kinematic properties of origami assemblages, realistic simulation models are needed. Although finite element (FE) simulations tend to provide a possible solution, the extended model building period and long computation time make them difficult to use at a preliminary design stage, where the instant feedback is desired. Several simplified models have been proposed to simulate the behavior of origami structures [710], and among them, two of the most widely used are rigid folding models and bar and hinge models.

The rigid folding models assume rigid panels and concentrate all deformation within creases. Numerical codes of this type are based on rigid folding constraints imposed by crease patterns. For example, one widely used necessary condition for rigid folding is that the rigid body rotations about a vertex sum to unity as R1R2, …, Rn = I, which was discovered by Kawasaki and Yoshida [11]. The detailed equations for expressing Ri can be found in Ref. [12]. Widely used numerical simulations of this type include Tachi’s rigid origami simulation code for triangulated general patterns [7] and quad origami patterns [8].

The bar and hinge models simulate an origami pattern using extensional bar elements and rotational springs so that they can further capture the deformation of panels. In a bar and hinge type model, the bar elements are used to capture in-plane behaviors such as panel stretching and shearing, while the rotational spring (hinge) elements are used to capture the out-of-plane behaviors such as panel bending and crease bending of origami structures. A more thorough introduction of this model is given in the Derivation of Stiffness Parameters section. Models of this type have been used to study small strain behaviors [13] as well as large strain behaviors of origami tessellations [9,10]. Moreover, they were also used in analyzing the mechanical characteristics of origami-inspired metamaterials [1,14]. A number of different variations of bar and hinge models were developed and studied in detail [9,15].

These two common origami simulation methods simplify creases as concentrated 1D rotational springs that have infinitesimal width (Fig. 1(a)). However, many active origami systems do not have these discrete hinges (such as those on our doors) because of limitations in the fabrication and assembling methods. For example, it is difficult to fabricate these discrete hinges due to the inherent 2D nature of microfabrication for small-scale origami [16]. These active origami systems are usually composed of softer compliant crease regions that can serve as actuators and panels that are stiffer and remain relatively rigid (Figs. 1(b) and 1(c)). In those applications, the existence of compliant creases introduces additional degrees of freedom and allows the structure to achieve richer deformations such as crease stretching and twisting.

Fig. 1
Fig. 1

The “smooth fold” model developed by Hernandez et al. [19,20] was one of the first to explicitly simulate the behavior of compliant creases without a FE discretization. This model is based on the rigid panel assumption and solves the deformed configuration by considering the geometric constraint imposed by the crease pattern. This model is capable of capturing accurate local geometry of compliant creases, including the curvature and relative extension, and saves $30%$ of the computation time when compared with FE simulations [19]. The model was used for designing a reconfigurable antenna [21] and a solar-powered underwater vehicle [22]. Although this model is a state-of-the-art in the modeling of compliant creases in origami assemblages and can capture the exact folded geometry, its derivation is complicated and it requires a long computation time. Another limitation of the model is that it cannot capture the flexibility of panels, because it uses the rigid panel assumption.

In this article, a relatively simple and computationally efficient bar and hinge model for compliant crease origami is developed and tested. By including additional bars and rotational springs in the crease area, the model is capable of capturing the additional torsional and extensional deformation in these compliant creases (see Figs. 2(a)2(c)). Because this model can capture the folding geometry of origami structures with sparsely located nodes, the model provides a computationally efficient solution for solving a large number of compliant crease origami structures within a short period of time at a preliminary design stage. Despite the model being relatively simple, it simulates the bistable behavior of a compliant crease origami well because the model captures the extensional and torsional deformations of the compliant creases. In origami vertexes that are nondevelopable such as a cone-like eggbox patterns or a hyperbolic-like arrangement of panels, bistability can only occur due to the extension of creases [2]. In other more traditional origami patterns, the inclusion of compliant creases gives additional degrees-of-freedom that can make the system potentially bistable or multistable. Beyond the ability to capture advanced stability behaviors, this compliant crease model can give more realistic geometrical representation of origami structures when simulating thickness and contacts between origami panels [23].

Fig. 2
Fig. 2

This paper first presents the formulation of the proposed compliant crease bar and hinge model, which is derived based on fundamental plate theories and the pseudo-rigid-body model. Next, a thorough verification is performed to compare the accuracy of the proposed model against FE simulations for both small strain stiffness and large strain stiffness. An early version of this compliant crease simulation was presented by the authors in a conference publication [24], while a finalized version of the formulation and access to the relevant computational codes are made available here. Moreover, this work presents simulation examples to show the efficiency and capability of the proposed model. Specifically these examples demonstrate that the proposed model can capture self-assembly and the bistability in compliant crease origami by capturing the creases bending, stretching, and twisting deformations.

## Derivation of Stiffness Parameters

In bar and hinge models, the bar refers to the extensional element that provides in-plane stiffness for shearing and stretching. The rotational spring (hinge) refers to the rotational element about a bar axis that provides out-of-plane stiffness for crease bending and panel bending (see Fig. 2(a)). We refer to Refs. [9,10] for details on how these elements are computationally formulated, and we focus on their stiffness definitions in this article. In this paper, the crease is different from the rotational spring, and it refers to the shaded area shown in Fig. 2(b), which is a collection of several bars, rotational springs, and nodes.

The meshing technique of the proposed model is developed based on the observed natural behavior of an origami assemblage with compliant creases. When a compliant crease is folded by applying pure moment on it, the curvature of the crease is evenly distributed across the width (Fig. 2(c)). The proposed model captures this even distribution of the curvature and the three-dimensional geometry by using three lines of rotational springs instead of just one. Two rotational springs connect the panels to the crease area, and two half-length rotational springs are located along the center of the crease region. When relative torsion occurs between two panels, the connecting crease deforms diagonally, as indicated in Fig. 2(d). This behavior is captured with the four additional rotational springs placed diagonally within the crease region. Notably, this meshing technique also captures the extensional stiffness of creases, which can play a significant role in analyzing bistable and multistable behaviors of origami [2].

### Derivation of Bar Areas.

The compliant crease bar and hinge model can be defined based on the following geometric parameters: the crease width W, crease length L, crease thickness tcrease, and panel thickness tpanel (Fig. 3(a)). Two indexes, crease aspect ratio and thickness ratio, are defined to quantify the geometrical characteristics of compliant creases. Crease aspect ratio is defined as: α = W/L, where W and L are crease width and crease length, respectively (see Fig. 3(a)). The thickness ratio is defined as: β = tpanel/tcrease, where tpanel and tcrease are the thickness of panels and creases, respectively. These two indexes describe geometric properties of compliant crease origami that affect the local and global stiffness of the system. Common origami designs with compliant creases typically have a crease aspect ratio α ≈ 0.1 and a thickness ratio β ≈ 3.0 [17,25].

Fig. 3
Fig. 3

Bar areas of the model are determined by matching the stiffness of the bar and hinge model with a theoretical plate model of the same bulk material. The bars within a crease region are divided into four different groups as shown in the Fig. 3(b). The area of A0 bars is later determined in the panel modeling section and is temporarily assumed to be infinite in this section. This assumption generally works well because panels are usually much thicker and stiffer than creases.

The area of A1 bars is calculated such that they produce the same stiffness as a plane-stress plate under a uniform tensile loading (Fig. 3(c)). The stiffness of the plate is calculated by applying a unit normal strain field $ϵ→=[100]T$, and the total extensional stiffness of the plate model is obtained as follows:
$ka,pl=Ecrease1−νcrease2tcreaseLW$
(1)
where the Ecrease and νcrease are the Young’s modulus and Poisson’s ratio of the crease region respectively. The axial stiffness of the compliant crease in the bar and hinge model is calculated with one additional assumption: the contribution of diagonal bars (A2 bars) is neglected because of their large inclination. This assumption tends to be valid because the crease length L is typically much larger than its width W. With this additional assumption, only four A1 bars are contributing to the axial stiffness, and thus, the stiffness of the bar and hinge model is calculated as follows:
$ka,BH=2EA1W$
(2)
By setting ka,BH = ka,pl, we solve the assigned area for A1 bars as follows:
$A1=Ltcrease2(1−νcrease2)$
(3)
Similarly, the area of diagonal bars (A2 bars) is set such that they produce the same shear stiffness as a plate with the same bulk material. A unit shear strain field $ϵ→=[001]T$ is applied, and the shear stiffness of the plate is calculated as follows:
$ks,pl=EcreasetcreaseL2W(1+νcrease)$
(4)
Because A0 bars are assumed to have infinite area, only the diagonal bars are contributing to the shear deformation in the bar and hinge model. Therefore, the shear stiffness of the bar and hinge model can be calculated as follows:
$ks,BH=2EcreaseA2(W2+L2)1.5L2$
(5)
Similarly, by setting ks,BH = ks,pl, we can solve the assigned areas for A2 bars. It is further assumed that A3 bars have the same area as the A2 bars. Overall, the A3 bars do not contribute significantly either to tensile or shear behaviors, because they are orthogonal to the edge of the crease region. Although this assumption is somewhat crude, it greatly simplifies the formulation and gives results that match FE simulations relatively well as will be demonstrated in the Model Verification sections. The areas of A2 and A3 bars are thus calculated as follows:
$A3=A2=(L2+W2)1.54(1+νcrease)tcreaseLW$
(6)

### Derivation of Rotational Spring Stiffness.

In this article, upper case Kspr denotes the length-normalized spring stiffness, and the lower case kspr is the total stiffness of the rotational spring, which is calculated as follows:
$kspr=KsprLspr$
(7)
where Lspr is the length of the rotational spring. A compliant crease within an origami structure can be treated as a compliant mechanism and modeled as a single rotational spring using the pseudo-rigid-body model described in Ref. [26] (Fig. 4(b)). The rotational stiffness of the corresponding pseudo-rigid-body crease model is given as follows:
$kcrease=EcreaseIW=EcreaseLtcrease312W$
(8)
In the proposed model, this single rotational spring is further divided into three lines of springs (the four Kspr1 in Figs. 4(a) and 4(b)), to better match the distributed curvature observed in real compliant creases. Because one spring is split into three lines of springs, the rotational stiffness of each spring should be increased by a factor of three such that it produces the same global stiffness as the one-spring pseudo-rigid-body model. Thus, the normalized rotational spring stiffness for these four Kspr1 springs is given as follows:
$Kspr1=3EcreaseIWL=Ecreasetcrease34W$
(9)
All aforementioned model parameters are derived analytically. However, it is difficult to solve the stiffness of the four diagonal rotational springs Kspr2 using the previous strategies, and thus, their stiffness is obtained by matching the results from FE simulations of compliant creases for a specific scenario. A factor f is used to scale the stiffness of Kspr2 with respect to Kspr1 as shown in Eq. (10). By matching the infinitesimal torsional stiffness of the bar and hinge model with the FE simulation for the α = 0.1 and β = 3.0 case, f = 4 is obtained and Kspr2 is calculated as follows:
$Kspr2=fKspr1=4×3EcreaseIWL=Ecreasetcrease3W$
(10)

It makes sense that this factor f should be greater than one, because the corresponding diagonal deformations have more concentrated curvatures as shown in Fig. 2(c), and thus, the diagonal springs should be stiffer. Details on the performance of this assigned f = 4 when other values of α and β are encountered will be demonstrated in the Model Verification section.

Fig. 4
Fig. 4

To sum up, the model parameters representing the crease stretching, shearing, and bending behaviors are derived analytically using the pseudo-rigid-body model and theoretical plate model. However, because of the simplicity of a bar and hinge representation, we were not able to derive the model parameters related with torsional behaviors of a crease analytically, and thus, an additional parameter f is used to match the torsional behavior of the proposed model to FE simulation result with a specific geometric setup with α = 0.1 and β = 3.0. These values are selected because they are representative of common origami structures [17,25]. This additional f factor only affects the stiffness of the rotational spring Kspr2 and thus does not influence the model behavior for crease stretching (controlled mainly by A1), crease shearing (controlled mainly by A2), and crease bending (controlled mainly by Kspr1). Generally speaking, we found that the suggested value f = 4 gives relatively good results over a large range of α and β values as will be demonstrated in Model Verification section.

### Modeling of Panels.

The area of the bars within panels is calculated using the generalized N5B8 model derived in Ref. [15] as shown in Eq. (11). The equation is formulated with the panel area S and total length of the bars Li within the panel and is derived such that the strain energy is preserved when the panel undergoes uniform in-plane dilation.
$Apanel=2Stpanel(1−ν)∑i=1Li$
(11)
The rotational spring stiffness within the panel is derived based on the knowledge that bending stiffness scales with the cubic of thickness and the material stiffness [27]. The cubic of the ratio between panel thickness and crease thickness and the ratio of Young’s modulus are used to scale the normalized crease stiffness to obtain the stiffness for panel springs as follows:
$Kpanel=Epaneltpanel3Ecreasetcrease3Kcrease$
(12)

More robust models for bar and hinge formulations of origami panels as those in Ref. [9] could be used, but they are derived based on fully restrained panel edges. For the purpose of this article, these two simple equations work relatively well because the deformation of panels is usually not the controlling factor that governs the global behavior of compliant crease origami as will be shown in the Simulation Examples section.

## Model Verification

In this section, the proposed compliant crease bar and hinge model is compared against FE simulations to verify its performance. This verification is divided into two parts: infinitesimal stiffness and large deformation stiffness. In Infinitesimal Stiffness Verification and Large Deformation Stiffness Verification sections, the elastic moduli of creases and panels are set to be Ecrease = Epanel = 2 GPa, while the Poisson’s ratios of both creases and panels are set as νcrease = νpanel = 0.3, which are realistic values of soft crease origami assemblages.

### Infinitesimal Stiffness Verification.

Four types of loading are checked in the infinitesimal stiffness verification. They are tensile, shear, bending, and torsional loading. Stiffness predictions from the compliant crease bar and hinge model are compared against FE simulations when varying α values or β values (defined in Model Derivation section). The value of α is changed by varying W while fixing L = 2.0 cm, and the value of β is changed by varying tpanel while fixing tcrease = 0.1 mm. These values were selected because they represent realistic small-scale origami structures.

The tensile verification for the bar and hinge model with respect to the FE simulation is shown in Fig. 5. Infinitesimal stiffness of both models is normalized by EcreaseLtcrease/W to give a dimensionless result. The proposed bar and hinge model is about $10−25%$ softer than the FE simulation for all explored geometric variations of the structure. This difference is attributed to the error in the panel region rather than the crease region as demonstrated in Fig. 5(e). The generalized N5B8 panel model as defined in Eq. (11) does not capture the Poisson’s effect correctly and underestimates the stiffness because it only matches the strain energy under uniform dilation. The total amount of deformation δ in the compliant crease region appears appropriate; however, the deformed shape does not match the FE simulation exactly.

Fig. 5
Fig. 5

The results of shear verification are summarized in Fig. 6. Similar to the tension case, the stiffness plotted is normalized by EcreaseLtcrease/W to give dimensionless results. Figures 6(c) and 6(d) show that the mismatch between the bar and hinge model and FE simulation is around $10−15%$, with the bar and hinge model predicting stiffer responses. This difference is most likely due to the error in the model formulation of panels rather than of creases (see Fig. 6(e)) due to the same reason discussed in the tensile verification.

Fig. 6
Fig. 6

Figures 7(a) and 7(b) illustrate the loading setup for infinitesimal bending verification, and the results are plotted in Figs. 7(c) and 7(d). The stiffness predictions from the two models match well for different α values when β ≈ 3.0 as illustrated in Fig. 7(c). When β becomes large (β > 5), the stiffness mismatch between the two models reaches an asymptotic value of about $25%$. This mismatch is due to the additional diagonal rotational springs Kspr2 that are introduced to capture torsional deformations of the crease. These diagonal springs soften the bending stiffness of the crease by introducing additional degrees-of-freedom. The softening is much less apparent at larger strains, as will be demonstrated in the Large Deformation Stiffness Verification section. When β is small, the mismatch is larger because bending of panels is significant in the global deformation, and it is not captured by the bar and hinge model.

Fig. 7
Fig. 7

The verification of torsional stiffness is summarized in Fig. 8. The results only match well when α ≈ 0.1 and differ significantly otherwise. When α is small, the contribution of torsional deformation from the panel is large, while when α is large, the pseudo-rigid-body model proposed in Eq. (8) does not hold. The bar and hinge model is not well suited for capturing these torsional deformations over a wide range of α values, and thus, we tune the factor f in Eq. (10) such that realistic results are obtained for a typical compliant crease region with α ≈ 0.1. The torsion results match well for different values of β when α = 0.1 as shown in Fig. 8(d).

Fig. 8
Fig. 8

Overall, the proposed bar and hinge model can predict infinitesimal stiffness of origami assemblages relatively well. The proposed model matches FE simulation well when α ≈ 0.1 and β ≈ 3.0, which are typical values found in real origami-inspired engineering structures. Exceeding these ranges tends to create errors in the simulations, but these errors are often still within acceptable ranges. Considering that the proposed compliant crease bar and hinge model only has about one tenth of the degrees-of-freedom when compared with a coarsely meshed FE simulation, this accuracy is acceptable.

### Large Deformation Stiffness Verification.

In the verification of large deformation stiffness, only pure bending and combined bending and torsion are studied because shear and tensile deformations are substantially smaller. Two methods are used to compare compliant crease bar and hinge model with the FE simulation: (1) relative error is used to measure the influence of β on the difference in response between the two models because force displacement curves generally overlaps for different values of β and (2) force displacement curves are used to demonstrate the influence of the α value, because relative error did not show clear correlation as α changes. The relative error is defined as follows:
$Error=1N∑i=1N(δiFEA−δiBH)2δiFEA$
(13)
where N is the total number of load increments, and δ*i is the predicted deformation of the ith load increment from the different models. The loading was set up such that the two models have the same force increment sequence. In the following verifications, α is varied by fixing L = 2.0 cm while changing W, and β is varied by fixing tcrease = 0.1 mm while changing tpanel.

Figure 9(a) illustrates the loading setup for verification of large deformation bending. Figure 9(b) shows that the bar and hinge model matches the FE model well when β is around 2–3. The error reaches an asymptotic value of $20%$ as β becomes large. This mismatch is also due to the additional diagonal springs Kspr2, which add additional flexibility as described in the Infinitesimal Stiffness Verification section. When β < 1.2, the results from the bar and hinge model become highly offset from the FE model. This offset is because the proposed bar and hinge model is not able to capture the bending behavior of the relatively thin panels. The large error when β < 1.2 is not of a concern because realistic origami-inspired structures with compliant creases usually have a larger β. Our verification indicates that the relative error is less affected by the α index and is around 5–10% for all α tested. Thus, the force displacement curves are plotted to better illustrate the influence of the α value. Figure 9(c) shows that the compliant crease bar and hinge model matches the FE simulation relatively well for all α values tested when β = 3.0.

Fig. 9
Fig. 9

Figure 10 shows the results for the verification of combined torsional and bending loading case. The trends found in the combined torsion and bending case are similar to those found in the pure bending case with large strain. When β > 2.0, the relative error between the FE simulation and the proposed bar and hinge model varies between $5%$ and $10%$. Larger errors are observed as the panels become very thin (β < 1.2). Similarly, the force displacement curves are plotted for various α values, and the relative errors for these curves are between $5%$ and $10%$.

Fig. 10
Fig. 10

## Simulation Examples

### Folding of a Pattern With Four Vertexes.

To demonstrate the capability and efficiency of the proposed model, we simulate the folding of a pattern with four vertexes by applying external loads on the origami structure as demonstrated in Ref. [19]. The equilibrium is traced with a modified general displacement controlled method used in Ref. [10]. Because the detailed dimensions of structure were not provided in the reference, the dimensions used in the bar and hinge model are assumed such that they provide a similar shape. The thickness of the panels is assumed to be five times the thickness of the creases, so that the structure behaves like a rigid panel origami. Although the bar and hinge model cannot recreate the exact smooth curvature at the folds, it provides the same global shape for the origami (Fig. 11(c)).

Fig. 11
Fig. 11

Simulation results and the computation time of the smooth fold model and FE simulation are obtained directly from Ref. [19]. The 20-s run-time for the bar and hinge model was recorded on a modern (2018) desktop computer with an i7-7700 processor. The proposed model was able to fold the structure with only about $5%$ of computation time needed for the smooth fold model. If only the global response and the geometry of the structure are of interest, using the proposed model can save a large amount of the computational effort. However, the proposed model does not have a high resolution of the local response, so the smooth fold model or the FE simulation would be more suitable if the local response is of interest.

### Self-Folding of a Flower Pattern.

Generally speaking, the folding speed of each crease within a rigid origami is different due to the geometrical constraints imposed by the pattern. However, real self-folding active origami usually actuate all creases simultaneously at similar rates and thus require special design of the pattern [16]. Here, we demonstrate how the proposed model can be used to approach this problem.

In this example, a flower-like pattern with 48 creases is used (see Fig. 12). This flower pattern has a diameter that is about 10 cm, and its creases can be self-folded to different angles and with different mountain or valley assignments. The panels have a Young’s modulus that is assumed to be 2000 MPa, and the creases have a modulus that is 20 MPa. The thickness of creases are 0.3 mm, and the thickness of panels are 1 mm. These dimensions are selected such that they represent realistic active origami structures. The self-folding is accomplished by changing the stress-free angle of each crease incrementally and tracing the equilibrium with a Newton-Raphson iteration method.

Fig. 12
Fig. 12

Both Figs. 12(a) and 12(b) attempt to fold the same structure into an umbrella-like configuration by assigning an alternating mountain–valley crease arrangement around the center vertex. These mountain–valley assignments give us a water bomb base. In the folding sequence presented in Fig. 12(a), the stress-free angle of all mountain folds is gradually changed to 0.5π, and the stress-free angle of all valley folds is gradually changed to π. This self-folding path meets the rigid folding requirements relatively well, and thus, the real folding angle follows the assigned stress-free angle closely. This folding path allows the structure to be successfully folded to the desired umbrella configuration with minimum strain energy developed as demonstrated in Fig. 12(d). However, if all creases are folded at the same rate, this pattern does not fold to the umbrella shape correctly (see Fig. 12(b)). In this folding path, the stress-free angles for both mountain folds and valley folds are gradually changed to 0.7π and thus have the same speed of folding. Because the folding path does not match well with the rigid folding constraints imposed by the pattern, the structure becomes trapped at an undesired configuration and the real folding angles diverge from the assigned stress-free angles extensively. This divergence generates a large amount of strain energy that is stored within the structure (Fig. 12(d)). In real practice, one can solve this problem by (i) adjusting the folding rate of each crease or (ii) assigning different folding sequences of each crease.

Figure 12(c) shows another possible folding path of the structure where each crease is folded such that the final geometry looks like a bowl. In this case, all active compliant creases are folded simultaneously, and their folding rates are the same. Because all creases are folded with the same rate, the folding path does not follow rigid folding requirements exactly. However, because the assigned folding path does not diverge from the rigid folding path significantly, the structure folds to its final configuration successfully despite building up strain energy in the system.

Although the pattern studied here has 48 creases and is a relatively complex problem, all simulations of the self-folding are finished within 1 min. This example further demonstrates the efficiency of the proposed model and shows that the model is suitable to study the self-folding process of compliant origami systems. The model allows researchers to investigate the possible design space rapidly and to learn about the global behaviors of the system before building more detailed FE models.

### A Bistable Fourfold Vertex.

In Fig. 13, we study the bistable behavior of a fourfold vertex presented in Ref. [2]. This four-fold origami vertex is a cone-like structure that experiences bistable snap-through behavior when loaded. Figure 13(a) shows the loading setup of the structure, where four nodes in the center are fixed in space while the loads are applied at the center of each panel. From the top view, these panels appear as 5 × 5 cm2 squares, and all compliant creases have the same 6 mm width. The cone shape is generated by moving the shared center nodes upward and the four outer free corners downward, which skews the square panel to form the missing angle as shown in Fig. 13(a). Thus, the actual panels are not square but are rhombus if flattened onto a plane. The thickness of the creases is 1 mm, and the thickness of the panels is 2.2 mm. These geometries are selected such that they represent the structure presented in Ref. [2]. The Young’s modulus of the crease material is set as 25 MPa, and the panel material is set as 2852 MPa. We directly use these values from Ref. [2] without altering them to recreate the test response recorded from Ref. [2] with our model. In principle, changing the material properties could eventually change the bistable behavior of these systems as well.

Fig. 13
Fig. 13

Figure 13(b) shows the force displacement curves recorded during the loading process, and Fig. 13(c) records the energy landscape. The force–displacement results obtained from the analyses are in good agreement with the experimental findings presented in Ref. [2], when comparing the initial peak force in the bistable snap-through. Furthermore, the plots show that the crease stretching energy dominates the energy barrier that forms the bistable snap-through behavior, again matching the observations presented in Ref. [2] well. Because conventional origami modeling techniques simplify the compliant creases as rotational hinges that have no width, they cannot capture the extensional deformations of the creases and thus cannot capture this type of bistable behavior. This example demonstrates that the proposed model provides a relatively fast and accurate simulation method to study bistable behaviors of origami where the bistability occurs due to the presence of compliant creases.

### Bistable Miura-ori Beams.

In this section, the bistable behavior of a Miura-ori beam with three units will be studied. The Miura-ori pattern is developable and can become bistable, for example, (1) if there is differential stiffnesses and rest angles in the fold lines [28], or (2) if applied constraints and loads cause reorientation of the geometry. Here, we explore the second case where a simple Miura-ori sheet is loaded to invert one of its vertices. To demonstrate this bistable behavior, we use a physical model with flexible polymer sheet creases and cardboard paper to serve as the panels (Fig. 14(a)). The geometry of the Miura-ori system is set as illustrated in Fig. 14(b), with a = b = 5 cm and $γ=80deg$. The thickness of the creases is 0.1 mm, and the thickness of the panels is 1 mm. The width of the creases is 6 mm for all compliant creases. We assume the panels have a Young’s modulus of 1 GPa, and the creases have a Young’s modulus of 2 GPa. When a pair of loads is applied to the outermost unit of the Miura cantilever (forces F in Fig. 14(a)), the outermost two panels bend over and reach another stable state. The left edge of the physical model is restrained to keep the angle $θ=20deg$ fixed throughout the simulation.

Fig. 14
Fig. 14

We perform numerical simulations for the same Miura beam folded to three different θ values. Figure 14(c) shows that by decreasing the folding angle θ we can obtain less obvious bistable behaviors, and at $θ=12deg$, we obtain a monostable structure. This behavior is expected because when the θ value is approaching zero, the Miura cantilever is not folded and is mostly flat. Thus, loading the structure using the pair of forces displayed in Fig. 14(a) is only bending the two creases elastically, which gives a monostable response. When the θ value becomes larger, the geometry of the folded Miura cantilever gives rise to the second stable configuration, where the last vertex has to be inverted into a new state that is not compatible with the initial rigid folding kinematics.

This is another example that demonstrates the benefits of having compliant creases modeled explicitly, especially when studying bistable behaviors in origami structures. Besides the previously mentioned stretching deformation of creases, this example demonstrates that relative torsional deformations between adjacent panels is also important for capturing the correct geometry of these bistable origami. Figure 14(d) shows how the proposed model captures this relative torsion between adjacent panels by including the additional diagonal Kspr2 springs (A2 bars) into the creases region.

## Concluding Remarks

This paper proposes an improved formulation for a bar and hinge model so that it can be used to study the behavior of origami assemblages with flexible compliant creases. Equations for assigning bar areas and spring stiffness are derived by matching the stiffness of the bar and hinge model with a theoretical plate model and the pseudo-rigid-body model for crease stretching, shearing, and bending behaviors. However, due to the simplicity of the proposed model, one-spring stiffness Kspr2 correlated to crease torsion cannot be obtained analytically, and thus, a factor f is used to scale it so that the torsional stiffness of our model matches the FE simulation for α = 0.1 and β = 3.0 case. We later confirm that this selected f = 4 factor performs relatively well for a much wider range of α and β values in Model Verification section.

The proposed model is checked against FE simulations, and the results match generally well under tensile, shear, bending, and torsion for infinitesimal deformation. The proposed model produces the best estimate of structural stiffness when the crease aspect ratio is α ≈ 0.1 and the panel thickness ratio is β ≈ 3.0. These ratios are representative of real origami systems with compliant creases. At large deformations, the behavior of the proposed model is not sensitive to the α value and is able to give closely matched predictions when β > 1.2.

Simulation examples were used to further demonstrate the capabilities and efficiency of the model. The model can capture the crease bending for three-dimensional folding simulations and can approximate self-assembly where folds are actuated. In self-assembly simulations, the model can predict if the desired shape will be reached with simultaneous actuation of all folds or if the successful folding requires separatedly controlled acutation of creases. Furthermore, the proposed model captures the contribution from torsional and extensional behaviors of creases, which are of great significance when studying the bistable or multistable behaviors of origami-inspired assemblages with nonnegligible crease width. Many origami patterns that are traditionally considered monostable are in fact multistable when the added flexibility of compliant hinges is taken into account. The proposed model can be used as a fast and relatively accurate modeling technique for preliminary exploration of origami-inspired assemblages with compliant creases. The geometry and the flexibility of compliant creases in origami structures provide a new basis for various engineering structures and metamaterials with special mechanical properties.

## Availability of Simulation Code

A matlab simulation code is developed to carry out the proposed simulation for origami with compliant creases. The code is open access and can be found at https://drsl.engin.umich.edu/ContactSimulation/. We are actively updating and debugging the simulation codes.

## Acknowledgment

The first author would like to acknowledge College of Engineering Dean’s Fellowship and College of Engineering Challenge Fellowship from College of Engineering, University of Michigan. The authors also acknowledge support from DARPA Young Faculty Award (Grant No. D18AP00071). This paper reflects the views and position of the authors and not necessarily those of the funding entities. The authors want to acknowledge the thoughtful comments and suggestions from all three reviewers. Their suggestions have improved the clarity and quality of this work substantially.

## References

1.
Filipov
,
E. T.
,
Tachi
,
T.
, and
Paulino
,
G. H.
,
2015
, “
Origami Tubes Assembled Into Stiff, Yet Reconfigurable Structures and Metamaterials
,”
,
112
(
40
), pp.
12321
12326
. 10.1073/pnas.1509465112
2.
Faber
,
J.
,
Arrieta
,
A. F.
, and
Studart
,
A. R.
,
2018
, “
Bioinspired Spring Origami
,”
Science
,
359
(
6382
), pp.
1386
1391
. 10.1126/science.aap7753
3.
Yasuda
,
H.
,
Tachi
,
T.
,
Lee
,
M.
, and
Yang
,
J.
,
2017
, “
Origami-Based Tunable Truss structures for Non-Volatile Mechanical Memory Operation
,”
Nat. Commun.
,
8
(
1
), p.
962
. 10.1038/s41467-017-00670-w
4.
Martinez
,
R. V.
,
Fish
,
C. R.
,
Chen
,
X.
, and
Whitesides
,
M.
,
2012
, “
Elastomeric Origami: Programmable Paper-Elastomer Composites as Pneumatic Actuators
,”
,
22
(
7
), pp.
1376
1384
5.
Lv
,
C.
,
Krishnaraju
,
D.
,
Konjevod
,
G.
,
Yu
,
H.
, and
Jiang
,
H.
,
2014
, “
Origami Based Mechanical Metamaterials
,”
Sci. Rep.
,
4
(
5979
), pp.
1
6
. 10.1038/srep05979
6.
Fang
,
H.
,
Chu
,
S.-C. A.
,
Xia
,
Y.
, and
Wang
,
K.-W.
,
2018
, “
Programmable Self-Locking Origami Mechanical Metamaterials
,”
,
30
(
15
), p.
1706311
7.
Tachi
,
T.
,
2009
, “
Simulation of Rigid Origami
,”
Origami 4
,
Lang
,
R. J.
, ed.,
,
Dec. 8–10, 2015
,
CRC Press
,
Boca Raton, FL
, pp.
175
187
.
8.
Tachi
,
T.
,
Domingo
,
A.
, and
Lazaro
,
C.
,
2009
, “
Generalization of Rigid Foldable Quadrilateral Mesh Origami
,”
Proceedings of IASS Symposium
,
,
Sept. 28–Oct. 2
, pp.
2287
2294
.
9.
Filipov
,
E. T.
,
Liu
,
K.
,
Tachi
,
T.
,
Schenk
,
M.
, and
Paulino
,
G. H.
,
2017
, “
Bar and Hinge Models for Scalable Analysis of Origami
,”
Int. J. Solids Struct.
,
124
(
1
), pp.
26
45
. 10.1016/j.ijsolstr.2017.05.028
10.
Liu
,
K.
, and
Paulino
,
G.
,
2017
, “
Nonlinear Mechanics of Non-Rigid Origami: An Efficient Computational Approach
,”
Proc. R. Soc. A
,
473
(
2206
), p.
20170348
. 10.1098/rspa.2017.0348
11.
Kawasaki
,
T.
, and
Yoshida
,
M.
,
1988
, “
Crystallographic Flat Origamis
,”
Mem. Faculty Sci., Kyushu Univ.
,
42
(
2
), pp.
153
157
. 10.2206/kyushumfs.42.153
12.
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
, pp.
458
460
.
13.
Wei
,
Z. Y.
,
Guo
,
Z. V.
,
Gudte
,
L.
,
Liang
,
H. Y.
, and
,
L.
,
2013
, “
Geometric Mechanics of Periodic Pleated Origami
,”
Phys. Rev. Lett.
,
110
(
21
), p.
215501
. 10.1103/PhysRevLett.110.215501
14.
Schenk
,
M.
, and
Guest
,
S. D.
,
2013
, “
Geometry of Miura-Folded Metamaterials
,”
PNAS
,
110
(
9
), pp.
3276
3281
. 10.1073/pnas.1217998110
15.
Liu
,
K.
, and
Paulino
,
G.
,
2018
, “
Highly Efficient Structural Analysis of Origami Assemblages Using the Merlin2 Software
,”
Origami 7, the 7th International Meeting on Origami in Science, Mathematics and Education (7OSME)
,
Oxford, UK
,
Sept. 5–7
, pp.
1167
1182
.
16.
Na
,
J.-H.
,
Evans
,
A. A.
,
Bae
,
J.
,
Chiappelli
,
M. C.
,
Santangelo
,
C. D.
,
Lang
,
R. J.
,
Hull
,
T. C.
, and
Hayward
,
R. C.
,
2015
, “
Programming Reversibly Self-Folding Origami With Micropatterned Photo-Crosslinkable Polymer Trilayers
,”
,
27
(
1
), pp.
79
85
17.
Bassik
,
N.
,
Stern
,
G. M.
, and
Gracias
,
D. H.
,
2009
, “
Microassembly Based on Hands Free Origami With Bidirectional Curvature
,”
Appl. Phys. Lett.
,
95
(
9
), p.
091901
. 10.1063/1.3212896
18.
Mao
,
Y.
,
Yu
,
K.
,
Isakov
,
M. S.
,
Wu
,
J.
,
Dunn
,
M. L.
, and
Qi
,
H. J.
,
2015
, “
Sequential Self-Folding Structures by 3D Printed Digital Shape Memory Polymers
,”
Sci. Rep.
,
5
, p.
13616
. 10.1038/srep13616
19.
Hernandez
,
E. P.
,
Hartl
,
D. J.
,
Akleman
,
E.
, and
Lagoudas
,
D. C.
,
2016
, “
Modeling and Analysis of Origami Structures With Smooth Folds
,”
Comput. Aided Des.
,
78
, pp.
93
106
20.
Hernandez
,
E. P.
,
Hartl
,
D. J.
, and
Lagoudas
,
D. C.
,
2017
, “
Kinematics of Orgiami Structures With Smooth Folds
,”
ASME J. Mech. Rob.
,
8
(
6
), p.
061019
.
21.
Hernandez
,
E. P.
,
Hartl
,
D. J.
, and
Lagoudas
,
D. C.
,
2017
, “
Analysis and Design of an Active Self-Folding Antenna
,”
Proceedings of the ASME 2017 International Design Engineering Technical Conference and Computers and Information in Engineering Conference
,
Cleveland, OH
,
Aug. 6–9
,
Paper No. DETC2017-67855
.
22.
Hur
,
D. Y.
,
Hernandez
,
E. P.
,
Galvan
,
E.
,
Hartl
,
D.
, and
Malak
,
R.
,
2017
, “
Design Optimization of Folding Solar Powered Autonomous Underwater Vehicles Using Origami Architecture
,”
Prooceedings of the ASME 2017 International Design Engineering Technical Conferences and Computers and Information in Engineering Conference
,
Cleveland, OH
,
Aug. 6–9
,
Paper No. DETC2017-67848
.
23.
Zhu
,
Y.
, and
Filipov
,
E.
,
2019
, “
An Efficient Numerical Approach for Simulating Contact in Origami Assemblages
,”
Proc. R. Soc. A
,
475
(
2230
), p.
20190366
. 10.1098/rspa.2019.0366
24.
Zhu
,
Y.
, and
Filipov
,
E.
,
2019
, “
Simulating Compliant Crease Origami With a Bar and Hinge Model
,”
IDETC/CIE
,
Anaheim, CA
,
Aug. 16–19
,
Paper No. DETC2019-97119
.
25.
Yoon
,
C.
,
Xiao
,
R.
,
Park
,
J.
,
Cha
,
J.
,
Nguyen
,
T. D.
, and
Gracias
,
D. H.
,
2014
, “
Functional Stimuli Responsive Hydrogel Devices by Self-Folding
,”
Smart Mater. Struct.
,
23
(
9
), p.
094008
.
26.
Howell
,
L.
,
2001
,
Compliant Mechanisms
,
John Wiley and Sons, Inc.
,
New York
.
27.
Timoshenko
,
S. P.
, and
Krieger
,
S. W.
,
1959
,
Theory of Plates and Shells
,
McGraw-Hill Education
,
New York
, p.
5
,
Eq. (3)
.
28.
Waitukaitis
,
S.
,
Menaut
,
R.
, and
van Hecke
,
M.
,
2015
, “
Origami Multistability: From Single Vertices to Metasheets
,”
Phys. Rev. Lett.
,
114
(
5
), p.
055503
. 10.1103/PhysRevLett.114.055503