A Bar and Hinge Model for Simulating Bistability in Origami Structures With Compliant Creases

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 [7–10], 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 R₁R₂, …, R_n = I, which was discovered by Kawasaki and Yoshida [11]. The detailed equations for expressing R_i 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.

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

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.

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

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 t_crease, and panel thickness t_panel (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: β = t_panel/t_crease, where t_panel and t_crease 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].

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 A₀ 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.

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.

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 K_spr2 and thus does not influence the model behavior for crease stretching (controlled mainly by A₁), crease shearing (controlled mainly by A₂), and crease bending (controlled mainly by K_spr1). 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.

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.

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 E_crease = E_panel = 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.

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 t_panel while fixing t_crease = 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 E_creaseLt_crease/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.

The results of shear verification are summarized in Fig. 6. Similar to the tension case, the stiffness plotted is normalized by E_creaseLt_crease/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.

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

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

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.

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 K_spr2, 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.

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%$⁠.

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

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.

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.

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.

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 cm² 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.

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.

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.

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 K_spr2 springs (A₂ bars) into the creases region.

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

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.

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.