Homogenization framework for rigid and non-rigid foldable origami metamaterials

TL;DR

This paper develops a homogenization framework modeling origami faces as plates to accurately predict complex deformations in non-rigid foldable origami metamaterials, surpassing existing rigid models in accuracy.

Contribution

It introduces a novel continuum homogenization approach using plate elements for origami, capturing non-rigid deformation modes and improving predictive accuracy.

Findings

Effective elastic properties depend nonlinearly on fold angle and crease stiffness.

The framework achieves up to 12.9% error compared to detailed simulations.

Existing models show errors up to 161%, highlighting the new method's superiority.

Abstract

Origami metamaterials typically consist of folded sheets with periodic patterns, conferring them with remarkable mechanical properties. In the context of Continuum Mechanics, the majority of existing predictive methods are mechanism analogs which favor rigid folding and panel bending. While effective in predicting primary deformation modes, existing methods fall short in capturing the full spectrum of deformation of non-rigid foldable origami, such as the emergence of curvature along straight creases, local strain at vertices and warpage in panels. To fully capture the entire deformation spectrum and enhance the accuracy of existing methods, this paper introduces a homogenization framework for origami metamaterials where the faces are modeled as plate elements. Both asymptotic and energy-based homogenization methods are formulated and implemented. As a representative crease pattern, we examine the Miura origami sheet homogenized as an equivalent Kirchhoff-Love plate. The results reveal that certain effective elastic properties are nonlinearly related to both the initial fold angle and the crease stiffness. When benchmarked with results from fully resolved simulations, our framework yields errors up to 12.9\%, while existing models, including the bar-and-hinge model and the rigid-panel model, show up to 161\% error. The differences in errors are associated with the complex modes of crease and panel deformation in non-rigid origami, unexplored by the existing models. This work demonstrates a precise and efficient continuum framework for origami metamaterials as an effective strategy for predicting their elastic properties, understanding their mechanics, and designing their functionalities.

Tables

Table 1. Table 1: Geometric and material parameters of Miura origami metamaterial.

Panel side length	$a$	mm	20
Panel side length	$b$	mm	20
Sector angle	$\gamma$	∘	60
Initial fold angle	${\theta }_0$	∘	30
Panel thickness	$t$	mm	0.13
Unit length crease stiffness	$K_{\text{cr}}/a$	N	0.1
Base material Young’s modulus	$E$	MPa	4000
Base material Poisson’s ratio	$\nu$		0.38

Table 2. Table F.1: Geometric and material parameters of the cardboard Miura specimen.

Panel side length	$a$	mm	20
Panel side length	$b$	mm	20
Sector angle	$\gamma$	∘	60
Initial fold angle	${\theta }_0$	∘	57
Panel thickness	$t$	mm	0.28
Unit length crease stiffness	$K_{\text{cr}}/a$	N	0.032
Base material Young’s modulus	$E$	MPa	7900
Base material Poisson’s ratio	$\nu$		0.3

Table 3. Table G.2: Geometric and material parameters of curve crease origami metamaterial.

Unit cell length	$L$	mm	36.06
Unit cell width	$W$	mm	33.28
Unit cell height	$H$	mm	8.66
Spline $x_2$-length	$V$	mm	11.09
Panel thickness	$t$	mm	0.13
Base material Young’s modulus	$E$	MPa	4000
Base material Poisson’s ratio	$\nu$		0.38

Equations

$$u_1(x_1,x_2,x_3)=u_{01}(x_1,x_2)-x_3{\partial }_{x_1}{u}_{03}(x_1,x_2)$$

$$u_2(x_1,x_2,x_3)=u_{02}(x_1,x_2)-x_3{\partial }_{x_2}{u}_{03}(x_1,x_2),$$

$$u_3(x_1,x_2,x_3)=u_{03}(x_1,x_2),$$

$$e_{11}={\partial }_{x_1}{u}_{01},$$

$${\kappa }_{11}=-{\partial }_{x_1}^2{u}_{03},$$

$$\left[\begin{array}{c}\mathbf{N}\\ \mathbf{M}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{A}}^H& {\mathbf{B}}^H\\ {\mathbf{B}}^H& {\mathbf{D}}^H\end{array}\right]\left[\begin{array}{c}\mathbf{e}\\ \mathbf{\kappa }\end{array}\right],$$

$${\bar{U}}^H=\frac{1}{2}{\left[\begin{array}{c}\mathbf{e}\\ \mathbf{\kappa }\end{array}\right]}^T\left[\begin{array}{cc}{\mathbf{A}}^H& {\mathbf{B}}^H\\ {\mathbf{B}}^H& {\mathbf{D}}^H\end{array}\right]\left[\begin{array}{c}\mathbf{e}\\ \mathbf{\kappa }\end{array}\right].$$

$$C_{ijk\ell }^{h\epsilon }(x_1,x_2,x_3)=\frac{1}{h^3}{C}_{ijk\ell }(\frac{x_1}{\epsilon },\frac{x_2}{\epsilon },\frac{x_3}{h}),$$

$${\partial }_{z_1}{\sigma }_{i1}^{\epsilon }+{\partial }_{z_2}{\sigma }_{i2}^{\epsilon }+\frac{1}{h}{\partial }_{z_3}{\sigma }_{i3}^{\epsilon }=0$$

$${\sigma }_{ij}^{\epsilon }=C_{ijk1}^{\epsilon }{\partial }_{z_1}{u}_k^{\epsilon }+C_{ijk2}^{\epsilon }{\partial }_{z_2}{u}_k^{\epsilon }+\frac{1}{h}{C}_{ijk3}^{\epsilon }{\partial }_{z_3}{u}_k^{\epsilon }$$

$$u_i^{\epsilon }=0$$

$${\sigma }_{1j}^{\epsilon }{n}_j=t_1/h,{\sigma }_{2j}^{\epsilon }{n}_j=t_2/h,{\sigma }_{3j}^{\epsilon }{n}_j=t_3$$

$$u_i^{\epsilon }(\mathbf{z})=u_i^{(0)}(\mathbf{x},\mathbf{y})+\epsilon u_i^{(1)}(\mathbf{x},\mathbf{y})+{\epsilon }^2{u}_i^{(2)}(\mathbf{x},\mathbf{y})+{\epsilon }^3{u}_i^{(3)}(\mathbf{x},\mathbf{y})+{\epsilon }^4{u}_i^{(4)}(\mathbf{x},\mathbf{y})+\dots$$

$${\sigma }_{ij}^{\epsilon }(\mathbf{z})={\epsilon }^{-4}{\sigma }_{ij}^{(-4)}(\mathbf{x},\mathbf{y})+{\epsilon }^{-3}{\sigma }_{ij}^{(-3)}(\mathbf{x},\mathbf{y})+{\epsilon }^{-2}{\sigma }_{ij}^{(-2)}(\mathbf{x},\mathbf{y})+{\epsilon }^{-1}{\sigma }_{ij}^{(-1)}(\mathbf{x},\mathbf{y})+\dots$$

$${\partial }_{y_j}\left[C_{ijpq}\right({\delta }_{pk}{\delta }_{q\delta }+\frac{1}{2}({\partial }_{y_q}{w}_p^{k\delta }+{\partial }_{y_p}{w}_q^{k\delta })\left)\right]=0\text{ on }Y,{\mathbf{w}}^{k\delta }\text{ is }(y_1,y_2)\text{-periodic,}$$

$${\partial }_{y_j}\left[C_{ijpq}\right(y_3{\delta }_{p\gamma }{\delta }_{q\delta }+\frac{1}{2}({\partial }_{y_q}{p}_p^{\gamma \delta }+{\partial }_{y_p}{p}_q^{\gamma \delta })\left)\right]=0\text{ on }Y,{\mathbf{p}}^{\gamma \delta }\text{ is }(y_1,y_2)\text{-periodic,}$$

$$A_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_Y{C}_{ijk\ell }({\delta }_{k\gamma }{\delta }_{\ell \delta }+\frac{1}{2}({\partial }_{y_{\ell }}{w}_k^{\gamma \delta }+{\partial }_{y_k}{w}_{\ell }^{\gamma \delta }))({\delta }_{i\alpha }{\delta }_{j\beta }+\frac{1}{2}({\partial }_{y_j}{w}_i^{\alpha \beta }+{\partial }_{y_i}{w}_j^{\alpha \beta })d\mathbf{y},$$

$$\mathbf{e}=\left[\begin{array}{cc}{e}_{11}& e_{12}\\ e_{12}& e_{22}\end{array}\right],2\mathbf{\mu }=\left[\begin{array}{cc}2{\partial }_{y_1}{w}_1& {\partial }_{y_2}{w}_1+{\partial }_{y_1}{w}_2\\ {\partial }_{y_2}{w}_1+{\partial }_{y_1}{w}_2& 2{\partial }_{y_2}{w}_2\end{array}\right],\mathbf{\kappa }=-\left[\begin{array}{cc}{\partial }_{y_1}^2{w}_3& {\partial }_{y_1}{\partial }_{y_2}{w}_3\\ {\partial }_{y_1}{\partial }_{y_2}{w}_3& {\partial }_{y_2}^2{w}_3\end{array}\right]$$

$$A_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_{Y^0}{C}_{\zeta \eta \theta \lambda }({\delta }_{\theta \gamma }{\delta }_{\lambda \delta }+\frac{1}{2}({\partial }_{y_{\lambda }}{w}_{\theta }^{\gamma \delta }+{\partial }_{y_{\theta }}{w}_{\lambda }^{\gamma \delta })+y_3{\partial }_{y_{\theta }}{\partial }_{y_{\lambda }}{w}_3^{\gamma \delta })$$

$$({\delta }_{\zeta \alpha }{\delta }_{\eta \beta }+\frac{1}{2}({\partial }_{y_{\eta }}{w}_{\zeta }^{\alpha \beta }+{\partial }_{y_{\zeta }}{w}_{\eta }^{\alpha \beta })+y_3{\partial }_{y_{\zeta }}{\partial }_{y_{\eta }}{w}_3^{\alpha \beta })d{y}_1^{\prime }d{y}_2^{\prime }.$$

$$A_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_{Y^0}\mathbf{C}({\mathbf{e}}^{0e,(\gamma \delta )}+({\mathbf{\mu }}^{\ast e,(\gamma \delta )}+y_3{\mathbf{\kappa }}^{\ast e,(\gamma \delta )})):({\mathbf{e}}^{0e,(\alpha \beta )}+({\mathbf{\mu }}^{\ast e,(\alpha \beta )}+y_3{\mathbf{\kappa }}^{\ast e,(\alpha \beta )}))d{\mathbf{y}}^{\prime },$$

$$A_{\alpha \beta \alpha \beta }^H=\frac{2}{L_1{L}_2}U.$$

$$A_{\alpha \beta \alpha \beta }^H+A_{\alpha \beta \gamma \delta }^H+A_{\gamma \delta \alpha \beta }^H+A_{\gamma \delta \gamma \delta }^H=\frac{2}{L_1{L}_2}U,$$

$$A_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}U-\frac{A_{\alpha \beta \alpha \beta }^H}{2}-\frac{A_{\gamma \delta \gamma \delta }^H}{2}.$$

$$D_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_Y{C}_{ijk\ell }(y_3{\delta }_{k\gamma }{\delta }_{\ell \delta }+\frac{1}{2}({\partial }_{y_{\ell }}{p}_k^{\gamma \delta }+{\partial }_{y_k}{p}_{\ell }^{\gamma \delta }))(y_3{\delta }_{i\alpha }{\delta }_{j\beta }+\frac{1}{2}({\partial }_{y_j}{p}_i^{\alpha \beta }+{\partial }_{y_i}{p}_j^{\alpha \beta }))d\mathbf{y}$$

$$D_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_{Y^0}\mathbf{C}({\mathbf{e}}^{0\kappa ,(\gamma \delta )}+({\mathbf{\mu }}^{\ast \kappa ,(\gamma \delta )}+y_3{\mathbf{\kappa }}^{\ast \kappa ,(\gamma \delta )})):({\mathbf{e}}^{0\kappa ,(\alpha \beta )}+({\mathbf{\mu }}^{\ast \kappa ,(\alpha \beta )}+y_3{\mathbf{\kappa }}^{\ast \kappa ,(\alpha \beta )}))d\mathbf{y},$$

$$D_{\alpha \beta \alpha \beta }^H=\frac{2}{L_1{L}_2}U.$$

$$D_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}U-\frac{D_{\alpha \beta \alpha \beta }^H}{2}-\frac{D_{\gamma \delta \gamma \delta }^H}{2}.$$

$$B_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_Y{C}_{ijk\ell }(y_3{\delta }_{k\gamma }{\delta }_{\ell \delta }+\frac{1}{2}({\partial }_{y_{\ell }}{p}_k^{\gamma \delta }+{\partial }_{y_k}{p}_{\ell }^{\gamma \delta }))({\delta }_{i\alpha }{\delta }_{j\beta }+\frac{1}{2}({\partial }_{y_j}{w}_i^{\alpha \beta }+{\partial }_{y_i}{w}_j^{\alpha \beta }))d\mathbf{y},$$

$$B_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}{\int }_{Y^0}\mathbf{C}({\mathbf{e}}^{0\kappa ,(\gamma \delta )}+({\mathbf{\mu }}^{\ast \kappa ,(\gamma \delta )}+y_3{\mathbf{\kappa }}^{\ast \kappa ,(\gamma \delta )})):({\mathbf{e}}^{0e,(\alpha \beta )}+({\mathbf{\mu }}^{\ast e,(\alpha \beta )}+y_3{\mathbf{\kappa }}^{\ast e,(\alpha \beta )}))d\mathbf{y}.$$

$$B_{\alpha \beta \gamma \delta }^H=\frac{1}{L_1{L}_2}U-\frac{A_{\alpha \beta \alpha \beta }^H}{2}-\frac{D_{\gamma \delta \gamma \delta }^H}{2}.$$

$$\{\begin{array}{rlrl}& {\partial }_{x_{\beta }}{N}_{\alpha \beta }+T_{\alpha }=0& & \text{in }{\Omega }^0,\\ & {\partial }_{x_{\alpha }}{\partial }_{x_{\beta }}{M}_{\alpha \beta }+{\partial }_{x_{\alpha }}{Q}_{\alpha }-T_3=0& & \text{in }{\Omega }^0,\\ & N_{\alpha \beta }=A_{\alpha \beta \gamma \delta }^H\frac{1}{2}({\partial }_{x_{\delta }}{\bar{u}}_{\gamma }+{\partial }_{x_{\gamma }}{\bar{u}}_{\delta })+B_{\alpha \beta \gamma \delta }^H(-{\partial }_{x_{\gamma }{x}_{\delta }}^2{u}_3^{(0)})& & \text{in }{\Omega }^0,\\ & M_{\alpha \beta }=B_{\gamma \delta \alpha \beta }^H\frac{1}{2}({\partial }_{x_{\delta }}{\bar{u}}_{\gamma }+{\partial }_{x_{\gamma }}{\bar{u}}_{\delta })+D_{\alpha \beta \gamma \delta }^H(-{\partial }_{x_{\gamma }{x}_{\delta }}^2{u}_3^{(0)})& & \text{in }{\Omega }^0,\\ & u_3^{(0)}=0,{\bar{u}}_{\alpha }=0& & \text{in }{\Gamma }^0,\end{array}$$

$$e_{11}=1,{\mathbf{\chi }}_{\text{node}}^0=\left[\begin{array}{c}{y}_1\\ 0\\ 0\end{array}\right];$$

$${\kappa }_{11}=1,{\mathbf{\chi }}_{\text{node}}^0=\frac{1}{2}\left[\begin{array}{c}2y_1{y}_3\\ 0\\ -y_1^2\end{array}\right];$$

Full text

\affiliation

organization=Department of Mechanical Engineering, McGill University, addressline=845 Sherbrooke Street West, city=Montreal, postcode=H3A0G4, state=Quebec, country=Canada

{highlights}

Homogenization platform for property prediction of origami metamaterials.

Asymptotic and energy-based homogenization with numerical implementation.

Equivalent continuum of origami metamaterials with high panel flexibility.

1 Introduction

Origami, the art of paper folding, enables the folding of a flat, thin sheet into a three-dimensional structure with remarkable kinematics and mechanical properties. Initially considered as an art, origami has so far inspired the development of multiple technologies for applications, such as soft robotic actuators (Chi et al., 2022; Li et al., 2019b; Yang et al., 2024; Xue et al., 2024), deployable structures (Dang et al., 2025; Jamalimehr et al., 2022; Overvelde et al., 2017; Melancon et al., 2021; Fang et al., 2024; Li et al., 2023), energy absorbers (Jamalimehr et al., 2026; Jiang et al., 2025; Almessabi et al., 2024), and foldable electronics (Zhang et al., 2022; Li et al., 2019a). Crease patterns with period tessellation enable the design of mechanical metamaterials, namely origami metamaterials, whose physical properties are primarily governed by the architecture of the pattern rather than the base material. Figure 1 shows some examples of origami metamaterials. A wide range of properties have been investigated thus far, such as energy absorption capacity (Gattas and You, 2015; Harris and McShane, 2021), wave propagation (Oudghiri-Idrissi and Guzina, 2022; Zhang and Rudykh, 2024; Yasuda et al., 2019), and tunable or re-programmable geometry (Mirzajanzadeh and Pasini, 2025; Nassar et al., 2022; Castle et al., 2016; Gao et al., 2023).

Despite the wealth of research in origami metamaterials, modeling their mechanical behavior still poses significant challenges. One of the earliest approaches, the rigid-panel model (Schenk and Guest, 2013; Wei et al., 2013; Hu and Liang, 2020), treats an origami metamaterial as an assembly of rigid panels connected with revolute joints forming a mechanism. The kinematics can thus be studied through the theory of mechanisms, and the structural mechanics through the force-displacement relation. Later, the bar-and-hinge model (Schenk and Guest, 2011; Filipov et al., 2017; Liu and Paulino, 2017; Overvelde et al., 2017) was proposed to include the role of panel deformation. The origami metamaterial was thus treated as a combination of a truss lattice and rotational springs (Lahiri and Pratapa, 2023). The former allows both extension in truss elements and revolution around pin joints, capturing the primary deformation modes at the global level. The bar-and-hinge model has been used to capture the overall shape change under large deformation (Xu et al., 2024, 2025). However, local stress and strain are often neglected or simplified by this method. For non-rigid foldable origami patterns, another popular option is to model their panels with plate Finite Elements (FE) (Hu et al., 2021; Filipov et al., 2015; Feng et al., 2022). This method gives an accurate description of the stress and strain distribution but is computationally demanding. The computational requirement becomes overwhelming when there is a large origami tessellation with a complex fold pattern. For this reason, the FE analysis of origami specimens is limited to a few unit cells of relatively simple patterns.

To bypass the fully resolved, detailed simulation of a periodic metamaterial, homogenization is often used to approximate the global response using an equivalent continuum. Multiple theories exist for homogenization, and are capable of simplifying the mechanical analysis of large, periodic tessellations, such as cellular origami metamaterials (Cheung et al., 2014; Lyu et al., 2021), origami sheets (Heimbs, 2013) and origami tubes (Turco et al., 2024), corresponding to 3D, 2D and 1D tessellations, respectively, for both the in-plane and out-of-plane mechanical properties. Xu et al. (2024, 2025) proposed an effective surface theory and its numerical implementation on nonlinear homogenization to characterize origami soft modes using the bar-and-hinge model. Cheung et al. (2014) applied periodic boundary conditions to an interleaved tube cellular structure to examine its elastic properties with classical dimensional scaling analysis. Vasudevan and Pratapa (2024) presented an energy-based homogenization framework with bar-and-hinge to determine the effective elastic constants of a couple-stress continuum representing origami sheets. Nassar et al. (2017) leveraged asymptotic homogenization to fit smooth surfaces with bar-and-hinge origami tessellations. Regardless of the method, the idea of homogenization proves to be a powerful tool in revealing the unique kinematics and mechanics of origami metamaterials.

Existing methods to homogenize the mechanical properties of origami metamaterials mainly capture rigid-foldable modes and rely on specific assumptions for folding. One is that creases and panels have contrasting properties, enabling the metamaterial to fold as a mechanism. Another is to treat the metamaterial as a non-reconfigurable structure that deforms in the panels rather than folding, similarly to a foldcore sandwich plate designed for load bearing and energy absorption (Zheng et al., 2022; Sturm et al., 2015; Zhang et al., 2021). Most works adopt either of these assumptions as they are easily implemented with the modeling methods reviewed above. However, origami systems may encompass a wide range of panel and crease characteristics. An example is soft origami actuators made of elastomer (Martinez et al., 2012; Jin et al., 2022). The soft and stretchable base material allows panels to bulge and wrinkle. Describing such mechanical responses is challenging for the rigid panel model and the bar-and-hinge model. At the same time, elastomeric creases may not show contrasting mechanical properties compared to the surrounding panels. Although this does not prevent the origami structure from folding into a wide range of motion, the assumption of rigid folding no longer applies. This is one among several other examples that call for the development of a property prediction method which is compatible with a broad variety of attributes in panels and creases so as to capture the mechanics of non-rigid foldable origami metamaterials. To address this gap, we investigate the previously uncharacterized deformation modes of origami metamaterials under varying crease stiffness.

This paper presents a homogenization framework for both non-rigid foldable and rigid foldable origami metamaterials that bridges the dichotomy between a mechanism and a structure. The framework predicts with accuracy when there is both folding kinematics and strain in the constituent material. As a demonstrative case, we examine a single-layered planar origami tessellation under stretching, shear, bending and twisting, which is of interest to a wide range of structural applications. Our method aims to homogenize origami metamaterials under the assumption that they can be treated as Kirchhoff-Love plates. We explore their linear elastic behaviors by computing the effective stiffness coefficients. To illustrate the versatility of the proposed framework, we apply two independent homogenization methods, namely asymptotic homogenization and energy-based homogenization. By comparing the results with those obtained with a detailed FE model of the corresponding origami tessellation, we determine the accuracy of effective properties for a wide range of the initial fold states. Next, we compare the accuracy of our method with the state-of-the-art, namely the bar-and-hinge model (Vasudevan and Pratapa, 2024) and the rigid-panel model (Wei et al., 2013). The homogenization methods can be straightforwardly extended to 3D for cellular origami metamaterials, such as the stacked Miura-origami and the Tachi-Miura polyhedron, but we present the case of a planar origami for brevity. Although the framework is general to a variety of origami patterns, we examine the Miura pattern as a representative pattern to demonstrate the procedure.

The remainder of the article is as follows: Section 2 presents the fundamentals for the homogenization methods, namely the asymptotic and the energy-based homogenization, which are suitable for both rigid and non-rigid origami patterns. The benchmark result of a Miura sheet is validated in Section 3. Section 4 investigates the roles of initial fold angle and crease stiffness and compares results with two well-known origami models in the literature. Section 5 highlights the differences between models in terms of crease and panel deformation under specific loads and boundary conditions. Section 6 concludes the study with final remarks.

2 Homogenization formulations for origami metamaterials

With a focus on a single origami sheet, we present the homogenization formulations for determining the effective elastic constants of an equivalent Kirchhoff-Love plate. A common approach of all homogenization methods is to take a representative region, which in our paper is a unit cell, and assume that its properties are equivalent to those of the entire metamaterial (Hollister and Kikuchi, 1992). This section presents two formulations, i.e., asymptotic and energy-based methods. They tackle the homogenization problem from dissimilar perspectives. The former approximates the field variables in the governing equation, such as the displacement, by an asymptotic expansion. The latter, on the other hand, approximates the in-situ stress or strain of the metamaterial with applied boundary conditions of a unit cell. Both methods are applicable to origami metamaterials under linear elastic deformation. The pattern should exhibit translational periodicity, with the unit cell size much smaller than the full tessellation. Finally, the unit cell height, length, and width should be on the same order of magnitude. By comparing both methods, we can appreciate their common and distinct traits, obtain the effective properties, and identify the strain energy density of the metamaterial under prescribed deformation modes.

2.1 Homogenized Kirchhoff-Love plate model for origami sheets

Figure 2 shows in blue the physical domain $\Omega$ of an origami sheet. Figure 2 shows its homogenized counterpart with equal length, width and height. The homogenization framework aims to find the effective properties of the homogeneous thin plate $\Omega$. Figure 2 introduces the 2D domain $\Omega^{0}$, which is both the neutral surface of the true origami sheet and the mid-surface of the homogenized one.

We select the neutral surface $\Omega^{0}$ such that it does not stretch or compress under pure bending (Vasudevan and Pratapa, 2024). Figure 2 presents the location of $\Omega^{0}$. The origami panels have a thickness of $t$. After folding, the sheet assumes an apparent global height of $h$, which increases with the folding angle $\theta_{0}$, shown in Figure 2. A Cartesian coordinate system is introduced with axes $x_{1}$, $x_{2}$, and $x_{3}$ where $x_{1}$ and $x_{2}$ lie on the neutral surface of the origami sheet and the mid-surface of the homogenized plate. The displacements of the homogenized plate are described by the displacements of its mid-surface as

[TABLE]

where $u_{01}$ and $u_{02}$ are the in-plane displacements of a point on the mid-surface, and $u_{03}$ is the out-of-plane deflection. Assuming small displacement gradients, we can express the strain components $\tilde{e}_{ij}$ in terms of the displacements $u_{i}$ as $2\tilde{e}_{ij}=\partial_{i}u_{j}+\partial_{j}u_{i}$. The following expressions adopt index notations with Greek letters $\alpha,\beta,\gamma,\delta=1,2$. The plate theory assumes that the strains $\tilde{e}_{13}$, $\tilde{e}_{23}$, and $\tilde{e}_{33}$ are zeros. The strains $\tilde{e}_{11}$, $\tilde{e}_{22}$, and $\tilde{e}_{12}$ can be reduced to functions of the mid-surface strains $e_{\alpha\beta}$ and curvatures $\kappa_{\alpha\beta}$ following $\tilde{e}_{\alpha\beta}=e_{\alpha\beta}+x_{3}\kappa_{\alpha\beta}$, where the mid-surface strains and curvatures read

[TABLE]

Adopting the Kirchhoff-Love plate model, the plate constitutive relations are (Reddy, 1996)

[TABLE]

where $\mathbf{N}$ is the resultant force vector, $\mathbf{M}$ is the resultant moment vector, $\mathbf{A}^{H}$ is the extensional stiffness coefficient matrix with components $A^{H}_{\alpha\beta\gamma\delta}$, $\mathbf{B}^{H}$ is the bending-extensional coupling stiffness coefficient matrix with components $B^{H}_{\alpha\beta\gamma\delta}$, and $\mathbf{D}^{H}$ is the bending stiffness coefficient matrix with components $D^{H}_{\alpha\beta\gamma\delta}$. The matrix of the stiffness coefficients, known as the ABD matrix, is symmetric with 21 independent components.

For homogeneous or laminated plates, there exist relationships between the stiffness coefficients $A^{H}_{\alpha\beta\gamma\delta}$, $B^{H}_{\alpha\beta\gamma\delta}$, $D^{H}_{\alpha\beta\gamma\delta}$ and the base material properties obtained from the Cauchy continuum theory (Reddy, 2006). However, this simplification does not apply to a nonhomogeneous medium such as an origami metamaterial with a general fold pattern. Rather, $A^{H}_{\alpha\beta\gamma\delta}$, $B^{H}_{\alpha\beta\gamma\delta}$, and $D^{H}_{\alpha\beta\gamma\delta}$ should all be computed independently (Vasudevan and Pratapa, 2024). The strain energy density of a homogeneous thin plate is

[TABLE]

The following introduces the formulations for asymptotic and energy-based homogenization, enabling us to determine the effective constants of a Kirchhoff-Love plate equivalent to the origami sheet under investigation.

2.2 Asymptotic homogenization

2.2.1 Theory

To apply the asymptotic homogenization theory to origami metamaterials, we assume that the apparent plate height $h$ is small and comparable to the length and width of a unit cell (Kalamkarov and Kolpakov, 1997). We start with a review of a modified asymptotic homogenization method for periodic thin plates (Caillerie and Nedelec, 1984; Kohn and Vogelius, 1984).

Consider the representative origami sheet shown in Figures 3 and 3. $\Omega_{h}$ denotes the domain of the sheet. A small parameter $h$ denotes the domain height of the equivalent homogenized plate. The upper, lower, and lateral surfaces are respectively $\Gamma_{h}^{+}$, $\Gamma_{h}^{-}$, $\Gamma_{h}^{\textup{lat}}$. $\Gamma_{h}^{+}$ and $\Gamma_{h}^{-}$ consist of all points $(x_{1},x_{2},x_{3})$ such that $(x_{1},x_{2})$ lies in the mid-surface $\Omega^{0}$. On $\Gamma_{h}^{+}$ $x_{3}=h$, while on $\Gamma_{h}^{-}$ $x_{3}=0$. $\Gamma_{h}^{\textup{lat}}$ has $(x_{1},x_{2})$ lying on the boundary of $\Omega^{0}$ and $x_{3}$ varying from 0 to $h$. There are $n$ unit cells along one given side of the sheet, with a total length of $\ell$. The origami pattern is periodic with a unit cell size of $\varepsilon=\ell/n$. Since the origami sheet is single-layered, its height $h$ and unit cell size $\varepsilon$ have the same order of magnitude. They are small parameters, meaning $h\to 0$ and $\varepsilon\to 0$ simultaneously as the origami tessellation becomes infinitely large.

To analyze the origami pattern, we rescale one of its unit cells $Y$ with a scaling factor $\varepsilon$. In addition to the coordinate system $\mathbf{x}=(x_{1},x_{2},x_{3})$ of the original problem, a new coordinate system $\mathbf{y}=(y_{1},y_{2},y_{3})$ is introduced, whose origin lies on the bottom left corner of the unit cell mid-surface. The macroscopic coordinates $\mathbf{x}$ for the origami sheet and the microscopic coordinates $\mathbf{y}$ for a unit cell are related through $y_{1}=x_{1}/\varepsilon$, $y_{2}=x_{2}/\varepsilon$, and $y_{3}=x_{3}/h$. The elastic stiffness of the material in the unit cell $C_{ijk\ell}$ is a function of $y_{1}$, $y_{2}$, and $y_{3}$. Before homogenization, the elastic stiffness of the full origami metamaterial $C_{ijk\ell}^{h\varepsilon}(x_{1},x_{2},x_{3})$ is periodic in the macroscopic coordinates $x_{1}$ and $x_{2}$ with its period identical to the size of the unit cell. It is assumed that

[TABLE]

where the superscripts $h,\varepsilon$ signify that the elastic stiffness is a function of the height and the periodicity. The indices $i,j,k,\ell$ range from 1 to 3. The following change of variables is now adopted to simplify the analysis: all functions of the origami sheet are made independent of $h$, i.e., $z_{1}=x_{1}$, $z_{2}=x_{2}$, and $z_{3}=x_{3}/h$. This results in a domain $\Omega$ for the origami sheet, an upper surface $\Gamma^{+}$, a lower surface $\Gamma^{-}$, and a lateral surface $\Gamma^{\textup{lat}}$. The elastic stiffness $C_{ijk\ell}^{h\varepsilon}(x_{1},x_{2},x_{3})$ becomes $C_{ijk\ell}^{\varepsilon}(z_{1},z_{2},z_{3})$, which is a function of the periodicity $\varepsilon$ but not a function of the height $h$. The full problem of elasticity of the origami sheet is

[TABLE]

where $\sigma_{ij}^{\varepsilon}$ is the stress tensor, $u_{i}^{\varepsilon}$ is the displacement, $n_{j}$ is the unit normal to the surface, and $t_{i}$ is the surface force. In Equations 6a and 6b, summation is only on the indices $i$, $j$, and $k$ but not on the fourth index. This is due to the $1/h$ term introduced by the change in variables from $x_{3}$ to $z_{3}$. We search for the solution to Equation 6 under the limit $\varepsilon\to 0$. The displacement is represented by the asymptotic expansion

[TABLE]

where each term $u_{i}^{(0)}$, $u_{i}^{(1)}$… are functions of both coordinates $\mathbf{x}$ and $\mathbf{y}$. These terms are smooth in $\mathbf{x}$ and periodic in $y_{1}$ and $y_{2}$. On the other hand, the stress tensor expands to

[TABLE]

where the terms $\sigma_{ij}^{(-4)}$, $\sigma_{ij}^{(-3)}$… have the same properties on $\mathbf{x}$ and $\mathbf{y}$ as the displacement terms. Applying the above expansions to Equation 6 gives an expression with various powers of $\varepsilon$. Grouping the terms with the same power of $\varepsilon$ yields a set of equations that simplifies to two microscopic problems.

The first one is

[TABLE]

where $\mathbf{w}^{k\delta}$ is a displacement field. The first term $\partial_{y_{j}}(C_{ijpq}\delta_{pk}\delta_{q\delta})$ represents the reaction force in the unit cell when it is extended or sheared to a strain equal to 1. $\mathbf{w}^{k\delta}$ is thus the solution of the microscopic problem corresponding to in-plane deformation modes.

The second microscopic problem is

[TABLE]

where $\mathbf{p}^{\gamma\delta}$ is a displacement field. The first term $\partial_{y_{j}}(C_{ijpq}y_{3}\delta_{p\gamma}\delta_{q\delta})$ represents the reaction force in the unit cell when it is bent or twisted to a curvature equal to 1. $\mathbf{p}^{\gamma\delta}$ is the solution of the microscopic problem corresponding to out-of-plane deformation modes.

Both microscopic problems are then rewritten in their weak forms before being solved through FEA, as detailed in the following section. The solutions $\mathbf{w}^{k\delta}$ and $\mathbf{p}^{\gamma\delta}$ are also known as the characteristic displacement fields (Hassani and Hinton, 1998).

After homogenization, the effective extensional stiffness coefficients are

[TABLE]

where $L_{1}$ and $L_{2}$ are the lengths of the unit cell along $y_{1}$ and $y_{2}$, respectively. Equation 11 applies to unit cells with a general shape. Since an origami unit cell consists of thin panels, we simplify the expression by assuming that the panels are Kirchhoff-Love plates. Their collective mid-surface $Y^{0}$ is shown in Figure 4. $Y^{0}$ is referred to the coordinates $(y_{1}^{\prime},y_{2}^{\prime})$. Note that the plates in a unit cell are not the homogenized plate introduced in Section 2.1. With a displacement field $\mathbf{w}$, the local strain field in the unit cell panels is $\mathbf{e}=\boldsymbol{\mu}+y_{3}\boldsymbol{\kappa}$ with

[TABLE]

where $\boldsymbol{\mu}$ is the in-plane strain and $\boldsymbol{\kappa}$, the curvature. The base material properties reduce to $C_{\zeta\eta\theta\lambda}$, with $\zeta,\eta,\theta,\lambda=1,2$, on the mid-surface $Y^{0}$ (Lewinski and Telega, 1988). Applying the above to Equation 11 gives

[TABLE]

Its shorthand expression is

[TABLE]

where $\mathbf{e}^{0e,(\alpha\beta)}$ denotes unit extension or shear strain in direction $\alpha\beta$, $\boldsymbol{\mu}^{*e,(\alpha\beta)}$ is the strain of the characteristic displacement $\mathbf{w}^{\alpha\beta}$ from the first microscopic problem in Equation 9, and $\boldsymbol{\kappa}^{*e,(\alpha\beta)}$ is the curvature of the characteristic displacement $\mathbf{w}^{\alpha\beta}$. The integral part of the expression is twice the strain energy of the unit cell. The relation between strain energy and effective properties applies to various types of structures, as seen later in our origami assembly of thin plates and torsional springs.

The strain energy is computed differently for diagonal terms $A^{H}_{1111},A^{H}_{2222},A^{H}_{1212}$ and off-diagonal terms $A^{H}_{1122},A^{H}_{1112},A^{H}_{2212}$. To find the diagonal terms $A^{H}_{\alpha\beta\alpha\beta}$, we sum the unit strain $\mathbf{e}^{0e,(\alpha\beta)}$ and the characteristic strain $\boldsymbol{\mu}^{*e,(\alpha\beta)}+y_{3}\boldsymbol{\kappa}^{*e,(\alpha\beta)}$ everywhere in the unit cell. The sum of strains has a corresponding strain energy $U$, which is expressed as a function of the effective stiffness coefficient, namely (Cheng et al., 2013)

[TABLE]

The rest of $A_{\alpha\beta\gamma\delta}^{H}$ is associated to two unit strains $\mathbf{e}^{0e,(\alpha\beta)}$ and $\mathbf{e}^{0e,(\gamma\delta)}$. Since our problem satisfies the principle of superposition, a linear combination of the two unit strains has a corresponding strain energy $U$, which gives

[TABLE]

where $A^{H}_{\alpha\beta\alpha\beta},A^{H}_{\gamma\delta\gamma\delta}$ are diagonal terms obtained from Equation 15, and $A^{H}_{\alpha\beta\gamma\delta}=A^{H}_{\gamma\delta\alpha\beta}$ from the symmetry of the stiffness matrix. The off-diagonal terms are expressed as

[TABLE]

Similarly,

[TABLE]

are the homogenized bending stiffness coefficients. Assuming the unit cell consists of thin plates, the expression simplifies to

[TABLE]

where $\mathbf{e}^{0\kappa,(\alpha\beta)}$ denotes unit curvature or twist in direction $\alpha\beta$, $\boldsymbol{\mu}^{*\kappa,(\alpha\beta)}$ is the strain of characteristic displacement $\mathbf{p}^{\alpha\beta}$ from the second microscopic problem in Equation 10, and $\boldsymbol{\kappa}^{*\kappa,(\alpha\beta)}$ is the curvature of the characteristic displacement $\mathbf{p}^{\alpha\beta}$. The unit cell strain energy $U$ is calculated from the sum of unit curvature $\mathbf{e}^{0\kappa,(\alpha\beta)}$ and the strain of corresponding characteristic displacement $\boldsymbol{\mu}^{*\kappa,(\alpha\beta)}+y_{3}\boldsymbol{\kappa}^{*\kappa,(\alpha\beta)}$. This strain energy relates to the diagonal terms as

[TABLE]

Note that Equation 20 differs from Equation 15 in the values of $U$. Similar to the in-plane case, there exists 3 linear combinations of unit curvature, each corresponding to a strain energy $U$. The off-diagonal terms are

[TABLE]

Another set of effective stiffness coefficients

[TABLE]

describes the coupling effect between in-plane and out-of-plane responses. It also simplifies to

[TABLE]

The integral part is twice the strain energy of a unit cell. All of $B_{\alpha\beta\gamma\delta}^{H}$ are off-diagonal and associated to the linear combinations of both in-plane strain and curvature. The sum of unit strains $\mathbf{e}^{0e,(\alpha\beta)}$, $\mathbf{e}^{0\kappa,(\gamma\delta)}$ and characteristic strains $\boldsymbol{\mu}^{*e,(\alpha\beta)}+y_{3}\boldsymbol{\kappa}^{*e,(\alpha\beta)}$, $\boldsymbol{\mu}^{*\kappa,(\gamma\delta)}+y_{3}\boldsymbol{\kappa}^{*\kappa,(\gamma\delta)}$ relates to the strain energy $U$. The coupling coefficients are

[TABLE]

With the above effective stiffness coefficients, we can now approximate the global response of an origami sheet as a thin plate. The governing equation Equation 6 homogenizes to

[TABLE]

where $\Omega^{0}$ is the mid-surface of the homogenized plate, $\Gamma^{0}$ is the boundary of the mid-surface, $N_{\alpha\beta}$ are the force resultants, $M_{\alpha\beta}$ are the moment resultants, $T_{1}$, $T_{2}$, and $T_{3}$ are averaged tractions defined as $T_{i}=\frac{h}{|Y|}\int_{\Gamma^{\pm}}t_{i}^{\pm}\mathop{\mathrm{d}\!}s$, $Q_{\alpha}$ are averaged bending moments defined as $Q_{\alpha}=\frac{h}{|Y|}\int_{\Gamma^{\pm}}y_{3}t_{\alpha}^{\pm}\mathop{\mathrm{d}\!}s$, $\bar{u}_{\alpha}$ are the homogenized displacements, and $u^{(0)}_{3}$ is the homogenized deflection.

2.2.2 Numerical implementation

Given the geometry of origami patterns, which can be complex, it is often not feasible to use closed-form expressions to perform asymptotic homogenization. Rather, a numerical approach becomes preferable. Our method here relies on FE analysis to find numerical solutions to the microscopic problems (Equations 9 and 10) and the effective stiffness coefficients (Equations 14, 19 and 23). The following procedure is fully implemented through ABAQUS python scripting, with python codes and input files available (Li, 2025). However, it is possible to generalize the procedure to any FE software package with minor modifications. Our implementation differs from the literature (Cai et al., 2014; Eskandari et al., 2024) in the strain energy calculation. Rather than calculating the strain energy from the nodal quantities during postprocessing, this step is performed fully using the selected FE software. This modification adds flexibility to the choice of element types. It is relevant to models of origami tessellations that require various types of plate elements, shell elements and torsional springs.

Taking the Miura sheet as an example, the procedure begins with identifying the geometry and material of the unit cell, specified in Table 1 along with Figure 5. We integrate this information into an FE model of the unit cell. We construct the unit cell (Filipov et al., 2015) with panels meshed as plate elements and creases as torsional springs. Since the choice of unit cell origin is arbitrary, in this example the bottom left corner of the Miura unit cell is selected as the center of a panel, shown in Figure 5. This choice simplifies the specification of boundary conditions in the microscopic problem.

Discretizing the microscopic problems exposed in Equations 9 and 10 leads to solving two systems of linear equations $\mathbf{K}\mathbf{w}^{\gamma\delta}=\mathbf{f}^{0}$ and $\mathbf{K}\mathbf{q}^{\gamma\delta}=\mathbf{f}^{0}$, where $\mathbf{K}$ is the stiffness matrix of the unit cell, $\mathbf{w}^{\gamma\delta}$ and $\mathbf{q}^{\gamma\delta}$ are the discretized characteristic displacement fields, and $\mathbf{f}^{0}$ is the initial strain loading. The latter has a physical meaning, which is nodal reactions induced by unit strains $e_{\gamma\delta}$ or unit curvatures $\kappa_{\gamma\delta}$. At the continuous level, the corresponding strain or curvature equals $1$ while all others are [math] everywhere in the unit cell. $\mathbf{f}^{0}$ is expressed as $\mathbf{f}^{0}=\int_{Y}\mathbf{B}^{T}\mathbf{E}\boldsymbol{\chi}^{0}\mathop{\mathrm{d}\!}\mathbf{y}$ (Cai et al., 2014), where $\mathbf{B}$ is the strain-displacement matrix, $\mathbf{E}$ is the constitutive matrix, and $\boldsymbol{\chi}^{0}$ is the displacement corresponding to $e_{\gamma\delta}$ or $\kappa_{\gamma\delta}$. $\boldsymbol{\chi}^{0}$ is obtained by constraining the displacement of each node, shown on the left of Figure 6. The constraints are expressed as

[TABLE]

where $\boldsymbol{\chi}^{0}_{\textup{node}}$ is the displacement vector of a node in the unit cell, and $(y_{1},y_{2},y_{3})$ are the coordinates of the node. Such constraints ensure that each displacement field $\boldsymbol{\chi}^{0}$ satisfies unit strain or curvature. The resulting unit cell with unit strain or unit curvature is presented in Figure 6.

In addition, the six deformation modes in Equation 26 have 15 linear combinations, as for instance

[TABLE]

The initial strain loading $\mathbf{f}^{0}$ is then found from

[TABLE]

Equation 26 and their linear combinations can be implemented by specifying a displacement boundary condition at each node. A subsequent static analysis in any FE software solves Equation 28. The initial strain loading $\mathbf{f}^{0}$ is extracted from the results as nodal reactions. The procedure is repeated 21 times ($6+15=21$) for the deformation modes in Equation 26 and their linear combinations.

The next step is to solve the microscopic problems

[TABLE]

where $\mathbf{f}^{0}$ is the initial strain loading obtained from the previous step, and $\boldsymbol{\chi}^{*}$ is the unknown characteristic displacement field. Different from the unit strains and curvatures, the characteristic displacements are assumed to be periodic in $y_{1}$ and $y_{2}$, as specified in Equations 9 and 10. This means that opposite unit cell boundaries have identical displacement. Because the Miura pattern tessellates in $(y_{1},y_{2})$, each unit cell has four boundaries, with the right boundary connected to the left boundary of an adjacent unit cell, and similar for the top and bottom boundaries. Accordingly the periodic boundary condition constrains the displacements of each pair of boundary nodes as

[TABLE]

where $\boldsymbol{\chi}^{*(\textup{l})}$, $\boldsymbol{\chi}^{*(\textup{r})}$, $\boldsymbol{\chi}^{*(\textup{t})}$, and $\boldsymbol{\chi}^{*(\textup{b})}$ are displacements on the left, right, top, and bottom boundaries, respectively.

The initial strain loadings $\mathbf{f}^{0}$ is applied to each node of the unit cell as a nodal force. Periodic boundary conditions are applied to the four sides of the unit cell in Figure 7 as nodal displacement constraints. Fixed boundary conditions to selected nodes are also enforced to prevent rigid body translation and rotation. Performing a static analysis in the FE software solves Equation 29. The procedure is repeated for all 21 initial strain loadings from the previous step to obtain the characteristic displacement fields $\boldsymbol{\chi}^{*}$. Figure 7 shows $\boldsymbol{\chi}^{*}$ corresponding to the deformation mode $e_{11}=1$.

The third step begins with finding the difference between the unit strain displacement field and the characteristic displacement field,

[TABLE]

$\tilde{\boldsymbol{\chi}}$ is applied to the unit cell to find its associated strain energy $U$. Figure 8 shows the displacement fields of six deformation modes as well as their strain energy. We specify $\tilde{\boldsymbol{\chi}}$ at each node of the unit cell through displacement boundary conditions. A static analysis is performed to find the total strain energy $U$. The analysis is repeated for the 6 deformation modes and their 15 linear combinations.

The fourth step is to find the effective stiffness coefficients $A_{\alpha\beta\gamma\delta}^{H}$, $B_{\alpha\beta\gamma\delta}^{H}$, $D_{\alpha\beta\gamma\delta}^{H}$ from Equations 15, 17, 20, 21 and 24.

The final step is the homogenized problem. The origami sheet is simplified into a thin plate. The homogenized plate $\Omega$, shown in Figure 2, is modeled with triangular or quadrilateral thin plate elements with a plate thickness of $h$. The effective elastic properties are specified in the FE model as elastic and anisotropic material properties. After applying appropriate loads and boundary conditions on the homogenized model, one last static analysis is performed, which solves the homogenized governing equation, see Equation 25.

A flowchart of the numerical implementation above is summarized in Figure D.4.

2.3 Energy-based homogenization

2.3.1 Theory

To apply the energy-based homogenization method, we assume that the origami metamaterial is represented by a unit cell, whose strain equals the average strain of the entire metamaterial. We apply periodic boundary conditions on the displacement to generate the desired strains, before solving the boundary value problem. Finally, we derive the effective stiffness coefficients from the unit cell strain energy.

A microscopic coordinate system $\mathbf{y}$ is introduced with its origin at the bottom left corner of the unit cell and $(y_{1},y_{2})$ axes on the mid-surface of the unit cell. Periodic boundary conditions are applied on the four sides of the unit cell. Specifically, the relative displacements are specified between opposite boundaries, generating specified strain or curvature in a given direction. The expressions for periodic boundary conditions, which generate in-plane strains, are (Wang et al., 2009)

[TABLE]

where $u_{\alpha}^{(\textup{r})}$, $u_{\alpha}^{(\textup{l})}$, $u_{\alpha}^{(\textup{t})}$, and $u_{\alpha}^{(\textup{b})}$ are displacements respectively on the right, left, top and bottom boundaries of the unit cell, $e_{\alpha\beta}$ are applied strains, and $L_{\beta}$ are the length of the unit cell along $y_{\beta}$. The periodic boundary conditions above differ from those of asymptotic homogenization (Equation 30) by a term $e_{\alpha\beta}L_{\beta}$, which generates an average strain of $e_{\alpha\beta}$.

The periodic boundary conditions corresponding to an average curvature $\kappa_{\alpha\beta}$ are (Vasudevan and Pratapa, 2024)

[TABLE]

where $x_{i}^{(\textup{r})}$, $x_{i}^{(\textup{l})}$, $x_{i}^{(\textup{t})}$, $x_{i}^{(\textup{b})}$ are positions respectively on the right, left, top and bottom boundaries of the unit cell.

The resulting strain energy density of a unit cell is

[TABLE]

where $Y$ is the mid-surface area, $\hat{\sigma}_{ij}$ is the local stress tensor, $\hat{e}_{ij}$ is the local strain tensor, and $y_{1},y_{2}$ are the local coordinates at the unit cell level. The energy density is evaluated per area of the unit cell mid-surface. The energy-based homogenization method requires that

[TABLE]

with $\bar{U}^{H}$ expressed in Equation 4. As per the symmetry of the stiffness coefficient matrix, 21 equations are required to compute the effective stiffness coefficients $A_{\alpha\beta\gamma\delta}^{H}$, $B_{\alpha\beta\gamma\delta}^{H}$, and $D_{\alpha\beta\gamma\delta}^{H}$. These equations are obtained from 21 unit strains, unit curvatures and their linear combinations applied to Equation 35. In practice, appropriate values are set for the strains $e_{\alpha\beta}$ and curvatures $\kappa_{\alpha\beta}$ to identify the corresponding deformation mode of the unit cell.

2.3.2 Numerical implementation

The numerical implementation of energy-based homogenization begins with constructing a model of the unit cell and identifying its boundaries. To obtain strains $e_{\alpha\beta}$, we apply constraints on each pair of nodes on opposite boundaries of the unit cell. For instance, to deform the unit cell into an average strain of $e_{11}=1$ on the boundaries of the unit cell, each node on the left boundary and its counterpart on the right boundary follow the constraints

[TABLE]

where $u_{i}^{(\textup{r})}$, $u_{i}^{(\textup{l})}$ are the right and left nodal displacements in $y_{i}$, and $y_{1}^{(\textup{r})}$, $y_{1}^{(\textup{l})}$ are the undeformed $y_{1}$ nodal coordinates. In parallel, each node pair on the top and bottom boundaries keeps its distance after deformation such that

[TABLE]

We apply the constraints in Equations 36 and 37 to the 3 translational DoFs of the corresponding node pairs. The constraints for in-plane deformation modes summarize as

[TABLE]

For out-of-plane deformation modes,

[TABLE]

where $y_{3}^{(\textup{r})}=y_{3}^{(\textup{l})}$ and $y_{3}^{(\textup{t})}=y_{3}^{(\textup{b})}$ due to geometric periodicity. In addition, we add appropriate boundary conditions to prevent rigid body motion. Figure D.4 compares the boundary conditions and loads of energy-based and asymptotic homogenization.

We numerically solve six boundary value problems with the above periodic boundary conditions enforced as nodal displacement constraints. Figure 9 shows the deformed unit cells. From the solutions of nodal displacements of each boundary value problem, we obtain the unit cell strain energy $U$. Figure 10 compares the total strain energy of each deformation mode.

In addition, we write the 15 linear combinations of the periodic boundary conditions in Equations 38, 39, 40, 41, 42 and 43. For example, the linear combination of $e_{11}$ and $e_{22}$ is known as biaxial extension. To find the corresponding boundary condition, we sum Equations 38 and 39 to obtain

[TABLE]

Similar to Equations 38, 39, 40, 41, 42 and 43, we apply Equation 44 as nodal displacement constraints. Solving the boundary value problem numerically gives the displacement field as well as the unit cell strain energy $U$ of the biaxial extension test. We then obtain the strain energy of the rest of the deformation modes.

The next step is to compute the effective stiffness coefficients $A_{ijk\ell}^{H}$, $B_{ijk\ell}^{H}$, and $D_{ijk\ell}^{H}$ from Equation 35 (Wang et al., 2009). Each strain energy $U$ computed from Equations 38, 39, 40, 41, 42 and 43 gives a diagonal term in Equation 4, while the linear combinations of the periodic boundary conditions give the off-diagonal terms.

Uniaxial extension, shear, bending, and twisting

Each set of periodic boundary conditions in Equations 38, 39, 40, 41, 42 and 43, corresponds to a total strain energy in the unit cell $U$. From the energy equivalence, we have

[TABLE]

where $L_{1}$ and $L_{2}$ are the lengths of the unit cell along $y_{1}$ and $y_{2}$, respectively.

Biaxial extension and other coupling terms

Each off-diagonal term of the effective stiffness matrix in Equation 4 corresponds to a strain or curvature in its row and another strain or curvature in its column. For example, $B^{H}_{1111}$ describes the coupling effect between uniaxial extension in $e_{11}$ and bending $\kappa_{11}$. Energy equivalence indicates that $B^{H}_{1111}$ is twice the strain energy density of the unit cell whose strains and curvatures are $e_{11}=\kappa_{11}=1$ and $e_{22}=e_{12}=\kappa_{22}=\kappa_{12}=0$. The strains and curvatures are imposed using a set of periodic boundary conditions, the linear combination of Equations 38 and 41. Then the strain energy $U$ is numerically computed. The effective properties are

[TABLE]

Note that the results from uniaxial extension, shear, bending, and twisting are required for $A^{H}_{\alpha\beta\alpha\beta}$ and $D^{H}_{\alpha\beta\alpha\beta}$ terms on the right-hand side of the equation above.

Finally, the homogenized stiffness matrix is applied to the homogenized model, as in the last step of Section 2.2.2.

A flowchart of the numerical implementation above is summarized in Figure D.5.

3 Validation

While an arbitrary origami sheet is described with 21 effective stiffness coefficients as in Equation 3, it is possible to reduce the number of constants by enforcing the symmetry of the origami pattern. For example, a Miura sheet exhibits centrosymmetry, meaning that a Miura sheet rotated by $180^{\circ}$ is identical to the original one. For centrosymmetric origami patterns, the following constitutive relation can be derived from the couple-stress theory (Vasudevan and Pratapa, 2024):

[TABLE]

where $\sigma_{\alpha\beta}$ are stresses, and $\mu_{\alpha\beta}$ are couple-stresses. The relationship involves only nine independent effective elastic constants, which are in-plane Poisson’s ratio $\nu_{12}$, Young’s moduli $E_{1}$ and $E_{2}$, shear modulus $G_{12}$, out-of-plane Poisson’s ratio $\nu^{b}_{12}$, bending moduli $M_{1}$ and $M_{2}$, and twisting modulus $T_{12}$. The out-of-plane Poisson’s ratio $\nu^{b}_{12}$ and $\nu^{b}_{21}$ are defined from the transverse curvatures, $\nu^{b}_{12}=-\kappa_{22}/\kappa_{11}$, and $\nu^{b}_{21}=-\kappa_{11}/\kappa_{22}$.

The effective elastic constants can be determined from the effective stiffness coefficients as

[TABLE]

where $\mathbf{S}$ is the compliance matrix in Equation 47, $h$ is the height of the origami sheet, and $[\mathbf{ABD}]$ is the homogenized matrix in Equation 3.

To set a baseline for comparison, we solved a detailed FE model of the origami sheet with identical parameters of the homogenized models, as listed in Table 1. Standard practice in homogenization applied to both periodic porous composites (Hollister and Kikuchi, 1992) and origami metamaterials (Vasudevan and Pratapa, 2024) shows that a detailed FE model is sufficient to validate the result from homogenization. To determine the elastic constants, we use a specimen with seven unit cells on each side of the sheet, adding up to 49 unit cells. On the other hand, for $G_{12}$ we examine a rectangular domain of 3 by 21 unit cells to avoid in-plane bending (Vasudevan and Pratapa, 2024). Although homogenization assumes an infinite tessellation of unit cells, a convergence test on the detailed FE model shows that 7 by 7 unit cells (or 3 by 21 for shear modulus) is a sufficient tessellation size. Simulations of six mechanical tests, shown in Figure 11, are required to determine all the elastic constants. To reduce the boundary effect, we measure the strain or curvature of the unit cell at the center of the tessellation and compute the reference elastic constants.

Figure 11 shows a uniaxial extension test with a roller boundary condition on the left boundary and a small displacement on the right boundary to generate an in-plane strain $e_{1}$. The resulting stress $\sigma_{1}$ and transverse strain $e_{2}$ are measured to obtain the Young’s modulus $E_{1}=\sigma_{1}/e_{1}$ and the Poisson’s ratio $\nu_{12}=-e_{2}/e_{1}$. Figure 11 shows the uniaxial extension along $x_{2}$, with a roller boundary condition on the bottom boundary and a displacement on the top. Measuring the stress $\sigma_{2}$ gives the Young’s modulus $E_{2}=\sigma_{2}/e_{2}$. Figure 11 is a shear test with the left boundary fully fixed and a displacement in $x_{2}$ applied on the right boundary. The shear modulus is $G_{12}=\sigma_{12}/(2e_{12})$. An additional roller boundary condition is required for the above tests to prevent global bending and twisting.

Figure 11 is a bending test with rollers on the left boundary and an applied bending moment $m_{2}$ evenly distributed on the right. The curvatures $\kappa_{1}$ and $\kappa_{2}$ are measured at the center of the specimen using the finite difference method. The bending modulus is $M_{1}=m_{2}/(hL_{2}\kappa_{1})$, where $H$ is the height, and $L_{2}$ is the width of the specimen in $x_{2}$-direction. The bending Poisson’s ratio is $\nu^{b}_{12}=-\kappa_{2}/\kappa_{1}$. Figure 11 is another bending test with rollers on the bottom boundary and a distributed moment $m_{1}$ on the top. We measure the curvature $\kappa_{2}$ at the center to find the bending modulus $M_{2}=m_{1}/(hL_{1}\kappa_{2})$. Figure 11 is a twisting test. Two corners on the diagonal are fixed in $x_{3}$, while the other two corners undergo a small displacement in $x_{3}$. The reaction force $F$ is measured at each corner, as well as the twist $\kappa_{12}$ at the specimen center. The twisting modulus is calculated as $T_{12}=F/(4h\kappa_{12})$. Additionally, the bending and twisting tests require boundary conditions to prevent rigid body motion.

Section B reports and discusses the errors with respect to the baseline results of the detailed FE models. Comparing the deviations of the plate FE asymptotic and energy-based homogenization, the bar-and-hinge model, and the rigid-panel model, we find that on average the errors are lowest for $\nu_{12}$, with the asymptotic plate FE \qty3.47%, the energy-based plate FE \qty3.20%, the bar-and-hinge \qty3.83%, and the rigid-panel \qty32.1%. The results of $G_{12}$ and $nu^{b}_{12}$ are similarly accurate. For $E_{1}$, $E_{2}$, $M_{1}$, and $M_{2}$, the averaged errors of the homogenized plate FE models are one third to one fifth those of other methods. In addition, the results for $T_{12}$ show that our framework is far more accurate (error of asymptotic plate FE \qty11.1%, energy-based plate FE \qty12.9%) than the bar-and-hinge model (error \qty161%).

To further validate our framework, we conducted a compression test on a Miura origami specimen to experimentally evaluate its effective Young’s modulus $E_{1}$. The specimen in Figure 12 is a 7 by 7 Miura tessellation made of cardboard. The fabrication procedure and specimen properties are detailed in Section F. The specimen is compressed with a universal testing machine between two flat plates with minimal friction. Figure 12 shows the experimental setup. We simulate the experiment with a homogenized model using energy-based method in Figure 12 and a detailed model in Figure 12. Both models have the same parameters as the specimen. The models share the same boundary conditions, with rollers on the bottom boundary and a uniform displacement of \qty10mm downwards on the top boundary. Comparing the results in Figure 12, we find the homogenized model in good agreement with both the detailed model and the experiment. In the linear elastic range, we compute the slope of the experimental curve at \qty1.98mm displacement, corresponding to a \qty1.04% strain. The error of the homogenized model is \qty34.7% compared with the experiment and \qty16.9% compared with the detailed model. The discrepancy between the homogenized model and the experiments can be attributed to manufacturing defects: the specimen is folded manually and its fold angle varies slightly from cell to cell. Other sources of error include a change of the crease properties due to repeated folding, the material anisotropy of the cardboard, and the friction during testing between the specimen and the holder.

4 Results

We can now apply both homogenization methods to determine the effective in-plane and out-of-plane properties of Miura origami sheets. Properties in the effective compliance matrix resemble those of a homogeneous thin plate made of an orthotropic material, with minor coupling between in-plane and out-of-plane responses. This observation confirms that the Miura pattern is centrosymmetric (Vasudevan and Pratapa, 2024).

To demonstrate the usage of the effective compliance matrices, we solve the global deformation of a Miura tessellation with seven unit cells on each side. Figure 13 shows a homogenized plate $\Omega$ with equal dimensions as the tessellation. The above compliance matrices are applied to the homogenized models as detailed in the last steps of Sections 2.2.2 and 2.3.2. The boundary condition on the left boundary is fixed in $u_{1}$ and $u_{2}$. There is a load of \qty0.4N evenly distributed on the right boundary, pointing [math] to the upper right. Figures 13, 13, 13 and 13 present the resulting displacement fields in the homogenized models. The detailed plate FE model, as the baseline, has a boundary condition of fixed $u_{1}$, $u_{2}$, and $u_{3}$ on the left boundary. There is a evenly distributed load of \qty0.4N on the right boundary, pointing [math] to the upper right. Figures 13 and 13 show the displacement fields. Both homogenized models agree well with the baseline in terms of the distribution of global displacement. There is discrepancy on the right boundary where the load is applied. Detailed model shows localized stretching of the panels that is not observed in the homogenized models.

The homogenized displacement fields are further validated against those from experiments. The cardboard Miura prototype and the simulations share the same geometric and material parameters, presented in Table F.1. The experimental setup ensures a fixed boundary condition on the left and a distributed load on the right. Figure 13 shows the deformed specimen. Its upper right corner undergoes a displacement magnitude of \qty10mm, which is identical to the detailed simulation at the same location. Section F details the experimental procedure.

To demonstrate the applicability range of the framework across fold patterns, we report in Section G the effective properties of a curve crease origami pattern adapted from the literature (Liu and James, 2024).

4.1 Role of initial fold angle

Since the architecture of a metamaterial governs its properties, we now examine how geometry tuning impacts the effective elastic constants of a Miura sheet. In this section, the selected geometric parameter is the initial fold angle $\theta_{0}$, as shown in Figure 14. It is measured as the dihedral angle between the panel and the $x_{1}$-$x_{2}$ plane. All four panels have the identical angle, a result stemming from the kinematics of Miura origami. $\theta_{0}$ represents the state of the pre-folded origami sheet, defined as the fold angle when the origami sheet is stress-free. We investigate the linear properties since the analysis does not consider large displacement or large strain. Among the geometric parameters listed in Table 1, we select $\theta_{0}$ as the varying parameter due to its strong influence on the folding kinematics. It also allows direct comparison with the bar-and-hinge model (Vasudevan and Pratapa, 2024) and the rigid-panel model (Wei et al., 2013), both of which have computed the effective properties of the Miura pattern under various fold angles. In this parametric study, $\theta_{0}$ varies from $1^{\circ}$ to $80^{\circ}$. Results from asymptotic and energy-based homogenization models are compared with those of the fully detailed FE model in Section 3 to determine result accuracy. After applying boundary conditions as shown in Figure 11, we characterize the detailed FE model in terms of its in-plane Poisson’s ratio $\nu_{12},$ Young’s moduli $E_{1}$, $E_{2},$ shear modulus $G_{12},$ out-of-plane Poisson’s ratio $\nu^{b}_{12},$ bending moduli $M_{1}$, $M_{2},$ and twisting modulus $T_{12}$. To reduce the boundary effect, we measure the above elastic constants using the unit cell at the center of the origami metamaterial.

Figures 15 and 15 compare the in-plane and out-of-plane kinematics of Miura origami through its Poisson’s ratio. The in-plane Poisson’s ratio is negative, a result that parallels the well-known observation of the auxetic response of the Miura sheet. A minimum occurs at $\theta_{0}=12^{\circ}$ with $\nu_{12}=-2.76$, while $\nu_{12}$ approaches zero when fully folded. The bending Poisson’s ratio is almost equal and opposite to the in-plane one, echoing the prediction of a pin-jointed truss framework (Schenk and Guest, 2011). These synchronized properties confirm that the Miura sheet behaves similarly to a single-DOF mechanism.

The in-plane moduli $E_{1}$, $E_{2}$, and $G_{12}$ share similar features at low to medium fold angles. For $\theta_{0}\approx 0^{\circ}$, the Miura sheet is almost fully flat, with properties reaching towards those of its base material, a thin sheet of Mylar. As the origami sheet folds, all effective moduli continuously decrease, because the relatively more flexible creases begin to play an increasingly influential role in the overall response. On the other hand, the kinematics of Miura pattern indicates that the total stiffness of the creases decreases for an increase in $\theta_{0}$ (Wei et al., 2013), a result that explains the trend observed on $E_{1},E_{2}$ in Figures 15 and 15. Figure 15 shows that the shear modulus $G_{12}$ is notably higher than the other two moduli. Previous work (Vasudevan and Pratapa, 2024) demonstrated that the shearing deformation mode involves mostly panel stretching, requiring a significant amount of strain energy. Finally, for $\theta_{0}>60^{\circ}$ there is a slight increase in $E_{2}$. This is explained by the alignment of panels along the $x_{2}$-direction at the highly folded state. In this case, loading along $x_{2}$ causes stretching in the panels and therefore an increase of $E_{2}$.

The out-of-plane moduli $M_{1}$, $M_{2}$, and $T_{12}$ follow a similar decreasing trend for $5^{\circ}<\theta_{0}<60^{\circ}$. However, at a small $\theta_{0}$ the bending modulus $M_{1}$ in Figure 15 shows a slight increase. This effect is similar to that of a corrugated board: the flutes makes it difficult to bend in one direction but easy in the other. The Miura pattern makes z-shaped flutes along $x_{1}$ that causes an increase in $M_{1}$, but as $\theta_{0}$ increases, this effect diminishes quickly since the folding kinematics has a strong influence on $M_{1}$. The ’flutes’ in the Miura sheet also makes $M_{1}$ higher than $M_{2}$ for most $\theta_{0}$. $M_{2}$ and $T_{12}$ show a decreasing to increasing trend as $\theta_{0}$ increases, as seen in Figures 15 and 15, which is qualitatively in agreement with the response of a truss-based lattice (Vasudevan and Pratapa, 2024).

4.2 Comparison with the literature

To illustrate the applicable range of the proposed framework, we compare results with other origami modeling methods in the literature, including the bar-and-hinge model (Vasudevan and Pratapa, 2024) and the rigid-panel model (Wei et al., 2013). Figure 16 summarizes all modeling methods under comparison.

As with the plate FE results presented in Section 4.1, the bar-and-hinge model has been homogenized with the energy-based method. The homogenized properties have been validated against a detailed bar-and-hinge model of 20 by 20 unit cells (with the exception of shear modulus, which requires a 4 by 50 tessellation). Both the bar-and-hinge models share the same parameters as the plate FE model, as listed in Table 1. The material properties are converted to the bar and spring properties as specified in previous work (Filipov et al., 2017). Another method in the literature is the rigid-panel model. It applies to the in-plane Poisson’s ratio $\nu_{12}$ and Young’s moduli $E_{1}$ and $E_{2}$. The shear modulus $G_{12}$ and all out-of-plane effective constants $\nu^{b}_{12}$, $M_{1}$, $M_{2}$, and $T_{12}$ are unbounded due to the assumption that the panels are perfectly rigid. The geometric parameters and crease stiffness are given in Table 1.

In Figure 17, the Poisson’s ratio $\nu_{12}$ closely aligns for $\theta>10^{\circ}$, with value increasing with $\theta$.

For $0^{\circ}<\theta<10^{\circ}$, the trend of both numerical models (bar-and-hinge and plate FE) differ from that of the rigid-panel model, with the former decreasing and the latter increasing. The monotonic increase in the rigid-panel model originates from the assumption that the origami metamaterial behaves as a 1-DoF mechanism. Its kinematic relation fully prescribes the mechanical response in Equation E.1. On the other hand, the two numerical models account for the base material properties. When $\theta$ approaches $0^{\circ}$, the origami metamaterial unfolds to a flat sheet of the base material, i.e., Mylar in this work. Correspondingly, the Poisson’s ratio from numerical models approaches that of Mylar. As $\theta$ varies, the base material properties and origami kinematics play a competing role in the metamaterial response. With the increase of $\theta$, the kinematic relationship gradually takes dominance. The Poisson’s ratio first decreases and then increases to align with the response of a pure mechanism, forming a local minimum at $\theta=10^{\circ}$ under the given parameters. Additionally, both the plate FE model and the bar-and-hinge model confirm that $\nu^{b}_{12}$ in Figure 17 is close to the negative of $\nu_{12}$ in Figure 17.

Figures 17 and 17 show that the in-plane moduli of all methods share the same trend. The Young’s moduli $E_{1},E_{2}$ of the plate FE model are lower than others. For example, at a given fold angle $\theta_{0}=30^{\circ}$, the moduli of the bar-and-hinge model are close to 1.5 times those of the plate FE. The discrepancy grows as $\theta_{0}$ becomes smaller and the Miura origami unfolds towards a flat sheet. The bar-and-hinge and the rigid-panel models align well with each other for most of the cases, when $10^{\circ}<\theta_{0}<80^{\circ}$. The shear modulus $G_{12}$ of plate FE and bar-and hinge models match in their values. The rigid-panel model is not included since it does not allow shearing.

The out-of-plane response involves bending and twisting. In Figure 17, the bending modulus $M_{1}$ of the plate FE model first increases then decreases as the Miura sheet unfolds, while the bar-and-hinge model has a monotonically decreasing trend. In Figure 17, $M_{2}$ of both models agree in their trends. Only when the sheet folds beyond $\theta_{0}>60^{\circ}$ their difference in values start to grow. $T_{12}$ of both methods also qualitatively agree in Figure 17. The result of the bar-and-hinge model is almost 3 times that of the plate FE, and this difference is consistent for most $\theta_{0}$.

4.3 Role of crease stiffness

This section investigates the role played by the constituent material in the effective elastic properties of an origami metamaterial. Some governing factors include the panel material properties and the crease formation process. Our paper adopts a simple crease model (Filipov et al., 2015) which neglects crease width and thickness and assumes a linear torsional stiffness evenly distributed along the crease line. The stiffness ranges from $0.1$ to \qty1.3N.mm for each \qty20mm long crease. To put numbers in context, the crease stiffness $K_{\text{cr}}$ of a typical origami metamaterial is smaller than or similar to the panel bending stiffness $K_{\text{p}}$. The panel bending stiffness is expressed as

[TABLE]

which is a function of the length of panel diagonal $D_{s}$, the panel thickness $t$, the base material Young’s modulus $E_{\text{base}}$ and Poisson’s ratio $\nu$, and the supplementary angle of a short diagonal corner $\alpha$ (Filipov et al., 2017). For the current panel geometry, $\sum\alpha=2\pi/3$. Models with varying crease stiffness are constructed, all with the parameters listed in Table 1.

Figure 18 presents the effective elastic constants for varying values of crease stiffness. Except for the in-plane and bending Poisson’s ratios, all constants are normalized by $E_{\text{base}}$. The crease stiffness is normalized by the panel bending stiffness as $K_{\text{cr}}/K_{\text{p}}$. A small ratio between crease and panel bending stiffnesses indicates that the origami sheet favors folding along crease lines, while a ratio close to 1 suggests that both crease folding and panel deformation have significant influence on the overall response of the origami sheet. Results of all relevant models, including the plate FE models with asymptotic and energy-based homogenization, the detailed plate FE model, the bar-and-hinge model, and the rigid-panel model, are all overlayed on the same plots.

Considering each effective elastic constants, Figures 18 and 18 shows that the crease stiffness has negligible influence on both in-plane and bending Poisson’s ratios $\nu_{12}$, $\nu^{b}_{12}$. This suggests that when the crease stiffness is smaller than or equal to the panel bending stiffness, the folding kinematics of Miura origami is consistent. In Figures 18 and 18, effective Young’s moduli $E_{1}$ and $E_{2}$ increase monotonically with the crease stiffness, demonstrating that the folding behavior of each crease has an influence on the global response of the metamaterial. Figure 18 shows that the crease stiffness does not affect the effective shear stiffness $G_{12}$ significantly, which confirms the observation that panel shearing is the main contributor to global shearing (Vasudevan and Pratapa, 2024). The bending moduli $M_{1}$ and $M_{2}$ increase with the crease stiffness as in Figures 18 and 18, meaning the stiffer the creases, the more difficult the origami sheet bends. The twisting modulus in Figure 18 has only a minor increase as the crease stiffness grows, if we consider the result of the detailed plate FE model. The creases do not influence the twisting modulus as much as the bending moduli, suggesting that panels contribute considerably to the global twisting behavior.

Although the global trends of the plate FE, the bar-and-hinge, and the rigid-panel models agree in general, there are some noticeable differences between their predictions. To begin with, the results of $\nu_{12}$, $\nu^{b}_{12}$, $G_{12}$, $M_{1}$, and $M_{2}$ show excellent agreement. Figure 18 shows that $G_{12}$ results of the bar-and-hinge model is closer to the baseline than the homogenized plate FE models. Discrepancy appears in $E_{1}$ and $E_{2}$ in Figures 18 and 18. The analytical results of the rigid-panel model overlap those of the bar-and-hinge model. They indicate that both $E_{1}$ and $E_{2}$ are linear functions of the crease stiffness. On the other hand, both the detailed and homogenized plate FE models predict that as the crease stiffness approaches the panel bending stiffness, results in the literature tend to overestimate $E_{1}$ and $E_{2}$, especially when the ratio $K_{\text{cr}}/K_{\text{p}}$ is above 0.4. Acknowledging that the rigid-panel model does not allow twisting, we compare the bar-and-hinge model with the plate FE models in Figure 18 and observe more deviation between them. The bar-and-hinge model indicates a linear relationship between the twisting modulus $T_{12}$ and the crease stiffness, with a slope larger than that of the plate FE models. This implies that to improve the accuracy of the effective elastic constants $E_{1}$, $E_{2}$, and $T_{12}$, we need to consider the relationship between crease and panel properties. Similarly, during mechanical simulations, it is also important to verify how significantly creases and panels influence the global behavior, while considering the applied loads and boundary conditions as another key factor.

5 Discussion

In this work, we present a homogenization framework for rigid and non-rigid origami metamaterials with panels modeled as plates. The framework applies to two independent methods, the asymptotic and the energy-based homogenization. The backbone of the framework unifying various homogenization methods is to apply appropriate boundary conditions to the unit cell, compute the resulting strain energy, and solve for the effective elastic constants using the constitutive relation. The majority of origami simulations in the literature include prescribed soft deformation modes according to the fold pattern. Regardless of the in-situ loads and boundary conditions, the mechanical response of the origami metamaterial always follows its soft modes. Our framework lifts this constraint by allowing the origami panels to deform as thin plates. We demonstrate the result accuracy over varying initial fold angles and crease stiffness.

To highlight the capabilities of our framework, we compare results with other existing origami models. Figures 17 and 18 show that a Miura sheet modeled with plate FE is more flexible than those with bar-and-hinge and rigid panels. The discrepancy between models is attributed to their applicable ranges of deformation modes. In general, the plate FE models (both detailed and homogenized) tend to be more compliant than the other models. This is true for both in-plane and out-of-plane responses. This result is attributed to a richer capture of deformation modes by the plate FE model, which allows more deformation modes in origami panels and creases than other models. As an example, Figure 19 shows that the rigid-panel model can only fold along its degree of freedom. The bar-and-hinge model, in addition to folding, can also stretch the panels or bend the panels along their diagonals. The plate FE model simulates all of the above and other non-rigid deformation modes such as bending and flexural behavior of creases and panels, along with localized twists at vertices. In contrast to the scaled deformation fields corresponding to unit strain or unit curvature in Figures 8 and 9, we present a realistic response of the detailed model, zoomed in on a unit cell in Figures 19 and 19. All figures show response in the linear elastic regime.

5.1 Stretching under constraints

To visualize the deformation modes of a Miura sheet, we apply a set of loads and boundary conditions in Figure 19 to the detailed plate FE and detailed bar-and-hinge models. There are forces applied at the top and bottom vertices, stretching the specimen in-plane. The left and right edges have roller boundary conditions. The bottom left corner is pinned to prevent rigid body motion. The applied boundary conditions prevent the origami sheet from unfolding along its degree of freedom, creating additional strain in the material. Zooming in on one of the unit cells at the center of the specimen, the plate FE model in Figure 19 shows a combination of bending and twisting in the panels that resembles warpage. Correspondingly, the initially straight creases deform into C-shapes that follow the curvature of the panels. At the vertex where corners of the panels meet, there is localized bending at the corners, which do not necessarily bend towards the same direction. Compared with the deformed bar-and-hinge unit cell in Figure 19 under the same loads and boundary conditions, the plate FE model exhibits a more significant change in shape, thus higher flexibility overall.

5.2 Twisting

Figure 19 examines the twisting deformation mode. The bottom left corner is pinned, and the top right corner is fixed in the out-of-plane direction. Loads are applied on the other two corners. The unit cell at the center deforms as in Figures 19 and 19. Twisting in each panel is captured with plate FE. In the bar-and-hinge model, the diagonal bar in each panel prevents twisting in certain directions. For example, the top-left panel in Figure 19 struggles to twist along the applied load due to the placement of its diagonal bar. Switching to another variation of the bar-and-hinge model with additional diagonal bars (Filipov et al., 2017) can potentially give more accurate results.

In summary, we compare the plate FE and the bar-and-hinge models under two sets of loads and boundary conditions to explain the discrepancies in the effective properties. Our result show consistent accuracy across a spectrum of initial fold angles and crease stiffnesses. We attribute the advantage of the proposed framework to a wide range of non-rigid deformation modes of the unit cell.

6 Conclusion

This work presents a homogenization framework to evaluate the linear elastic properties of an origami metamaterial through an effective computational model. The framework can model both non-rigid and rigid-foldable origami. It presents and compares two homogenization methods, asymptotic and energy-based homogenization. For demonstration, we implement both on a non-rigid Miura origami sheet and find its equivalent continuum model as a Kirchhoff-Love plate. A step-by-step procedure is provided to implement both homogenization methods. The results are benchmarked with a fully resolved numerical model of a finite Miura origami tessellation. The role of the initial fold angle is investigated through a parametric study. Both the base material properties and the folding kinematics have an influence on the effective properties. The variation of effective properties is also compared with those obtained from the rigid-panel and bar-and-hinge methods. Averaging the errors across a wide range of fold angles, the effective shear modulus has errors within \qty8% for all models. The effective twisting modulus shows considerable difference, with the error of the bar-and-hinge model 12 times that of the homogenization framework presented here. An origami sheet modeled with plates is more flexible than the literature, a result attributed to the variety of non-rigid deformation modes that plate elements can capture. Finally, a parametric study on crease stiffness demonstrates that both the literature and the proposed framework are sufficiently accurate for the effective in-plane and out-of-plane Poisson’s ratios, the shear modulus, and the bending moduli, whereas the proposed framework is recommended for the effective Young’s moduli and the twisting modulus.

Although our paper examines the basic Miura origami folded from a sheet of isotropic material, our framework applies to a wide range of origami patterns and base materials. Additionally, it can be easily implemented with the aid of any FE software package. The framework can be used for various applications, including origami-based soft robots, foldable antennas, and aircraft structures with periodicity. It is possible to extend the framework to cellular origami metamaterials such as the stacked Miura origami. In the context of Continuum Mechanics, it is straightforward to find the 3D analogy to our in-plane effective stiffness coefficients $A^{H}_{\alpha\beta\gamma\delta}$ following Section 2.2 or Section 2.3. Future work may involve extending the framework to nonlinear homogenization to incorporate the analysis of the curved-crease origami, such as the Miura lens-box pattern  (Mirzajanzadeh and Pasini, 2025).

7 Methods

This manuscript is accompanied with a set of ABAQUS and Python scripts recovering the central calculations in this work (Li, 2025).

8 Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

9 Acknowledgments

XL acknowledges support from Fonds de Recherche du Québec - Nature et technologies (Grant No. 348668). DP acknowledges the Canada Research Chair Program (Grant No. 258679), and the Natural Sciences and Engineering Research Council of Canada (Grant No. 208241). ML acknowledges the Natural Sciences and Engineering Research Council of Canada (Grant No. 401559).

XL gratefully acknowledges Morad Mirzajanzadeh and Filippo Agnelli for helpful discussions on the methodology and results.

A Mesh convergence test

To verify the proposed homogenization framework, we compute the effective elastic constants of a single-layered Miura origami sheet. Mesh convergence is checked for both asymptotic and energy-based homogenization. The results are first compared with a fully detailed FE model consisting of 49 unit cells. Next, we consider the limiting case of very stiff origami panels and compare it with an analytical model from the literature (Wei et al., 2013).

The geometric and material parameters of the origami metamaterial is given in Table 1. Mesh convergence of asymptotic and energy-based homogenization methods are compared in Figure A.1 at a crease stiffness of \qty0.05N. We increase the number of elements along each side of the rhombic panel and plot the variation of $x_{1}$-direction Young’s modulus $E_{1}$. The mesh convergence plot shows that the energy-based homogenization method converges at 16 elements on each side of the panel, while asymptotic homogenization requires 28 elements. In subsequent analysis, we ensure that both homogenization methods are converged with a mesh density of 28 elements per side.

B Error analysis

We consider the detailed plate FE model as the baseline. We compare its effective elastic constants with those of two homogenized plate FE models, a bar-and-hinge model (Vasudevan and Pratapa, 2024), and a rigid-panel model (Wei et al., 2013). In addition to the varying parameter of fold angle, which is in the range of $(0^{\circ},80^{\circ})$, all models share identical geometry and crease stiffness, as summarized in the first six columns of Table 1. The rigid-panel model differs from the rest in terms of panel flexibility. Due to the nature of the rigid panels, the in-plane shear modulus and all out-of-plane properties are unbounded. Only the in-plane Poisson’s ratio and Young’s moduli are compared with those of other models. The relative error between each model and the baseline is quantified as

[TABLE]

where $P$ represents any of the elastic constants $\nu_{12}$, $E_{1}$, $E_{2}$, $\nu^{b}_{12}$, $M_{1}$, $M_{2}$, and $T_{12}$. Their errors under varying initial fold angles $\theta_{0}$ are presented in Figure B.2.

The error over the fold angles is averaged as

[TABLE]

where $\theta_{0}^{\min}=0^{\circ}$ and $\theta_{0}^{\max}=80^{\circ}$. The maximum error in the range of fold angles between $0^{\circ}$ and $80^{\circ}$ is defined as

[TABLE]

Errors of all elastic constants are compared in Figure B.3.

On average, Poisson’s ratios $\nu_{12}$, $\nu^{b}_{12}$ and shear modulus $G_{12}$ are consistently accurate across all models. For the rest of the properties, the asymptotic and energy-based plate FE methods have relatively consistent ranges of error. They show an advantage in terms of accuracy when calculating the twisting modulus $T_{12}$ in particular. The maximum errors of asymptotic homogenization and energy-based homogenization are both within \qty53%.

C Effective compliance matrices of Miura origami

To provide an intuitive description of the mechanical response of a Miura sheet, we report the effective compliance matrix, which is the inverse of the stiffness coefficient matrix in Equation 3. Small numerical values in the matrix have been simplified to 0 because they are 2 to 4 orders of magnitude smaller than other entries, reflecting the planes of symmetry and hence anisotropy of the unit cells. The compliance matrix from asymptotic homogenization reads

[TABLE]

while for energy-based homogenization,

[TABLE]

The models of both matrices use geometric and material parameters in Table 1. The units in the matrix are \unit[]MPa^-1 for in-plane terms, \unit[]mm/N for in-plane and out-of-plane coupling terms, and \unit[]N^-1 for out-of-plane terms.

D Flowcharts of homogenization framework

D.1 Asymptotic homogenization workflow

Inputs

Unit cell geometry and material properties.

Step 1

Prescribe deformation fields $\boldsymbol{\chi}^{0}$ corresponding to unit strain or curvature according to Equations 26 and 27 on every node of the unit cell. Compute the initial strain loading $\mathbf{f}^{0}$.

Step 2

Apply periodic boundary conditions on the unit cell boundaries according to Equation 30. Numerically solve the microscopic problem to obtain the characteristic displacement fields $\boldsymbol{\chi}^{*}$.

Step 3

Find the difference between the unit strain displacement field $\boldsymbol{\chi}^{0}$ and the characteristic displacement field $\boldsymbol{\chi}^{*}$. Apply the difference as nodal displacements on the unit cell and compute its strain energy $U$.

Step 4

Compute the effective stiffness coefficients $A_{\alpha\beta\gamma\delta}^{H}$, $B_{\alpha\beta\gamma\delta}^{H}$, $D_{\alpha\beta\gamma\delta}^{H}$ from the unit cell strain energy $U$.

Homogenized problem

The full origami metamaterial is simplified into a thin plate with anisotropic material properties $A_{\alpha\beta\gamma\delta}^{H}$, $B_{\alpha\beta\gamma\delta}^{H}$, $D_{\alpha\beta\gamma\delta}^{H}$.

D.2 Energy-based homogenization workflow

Inputs

Unit cell geometry and material properties.

Step 1

Apply periodic boundary condition to the unit cell boundaries according to Equations 38, 39, 40, 41, 42 and 43.

Step 2

Numerically find the unit cell strain energy $U$ of each deformation mode in Step 1.

Step 3

Compute the effective stiffness coefficients $A_{\alpha\beta\gamma\delta}^{H}$, $B_{\alpha\beta\gamma\delta}^{H}$, $D_{\alpha\beta\gamma\delta}^{H}$ from the unit cell strain energy $U$.

Homogenized problem

The full origami metamaterial is simplified into a thin plate with anisotropic material properties $A_{\alpha\beta\gamma\delta}^{H}$, $B_{\alpha\beta\gamma\delta}^{H}$, $D_{\alpha\beta\gamma\delta}^{H}$.

Figures D.4 and D.5 present a summary and comparison of the numerical implementation of asymptotic homogenization and energy-based homogenization, including the required loads and boundary conditions, examples of deformation modes and examples of unit cell strain energy.

E Simulation methods for detailed origami metamaterials and unit cells

This section characterizes the linear elastic properties of an origami metamaterial, taking the Miura pattern as an example. The pattern tessellates in two directions, creating a planar metamaterial in the form of a single-layered sheet. Figure 2 shows the geometric parameters that defines the origami fold pattern and the initial fold state. $a$ and $b$ denote the panel side lengths, $\gamma$ is the sector angle, $\theta$ is the fold angle, and $t$ is the panel thickness. The material properties in a unit cell is characterized separately for the creases and the panels. The creases are simplified into torsional springs with a torsional stiffness of $K_{\textup{cr}}$ per crease. The panels are made of a linear elastic material with Young’s modulus $E$ and Poisson’s ratio $\nu$. We select their values as those of Mylar, a material widely used to fabricate origami specimens in the literature (Filipov et al., 2017). The geometric and material properties are summarized in Table 1.

Some of the most widely-adopted origami simulation methods include the rigid-panel, the bar-and-hinge, and the plate FE models. Sections E.1, E.2 and E.3 present each of their formulations.

E.1 Rigid-panel model

Assuming that the kinematic deformation modes are dominant in the structural behaviors of the metamaterial, the simplest rigid-panel model requires that the origami panels acts as rigid bodies, and creases as frictionless hinges. When the origami metamaterial folds, its panels rotate around each other following the creases. The Miura origami resembles a 1 degree-of-freedom (DoF) mechanism. Since the panels are perfectly rigid, it is only possible to stow and deploy in-plane. Taking the fold angle $\theta$ as the only DoF, we describe the in-plane Poisson’s ratio as (Schenk and Guest, 2013):

[TABLE]

Elastic responses appear when we consider the fold stiffness. By assuming that all elastic energy is stored in the creases, we simplify the folding process to rigid panels rotating around elastic hinges. The effective elastic constants in the two in-plane directions are (Wei et al., 2013):

[TABLE]

where $\rho$, chosen as the DoF in (Wei et al., 2013), is the dihedral angle between adjacent panels. Due to panel rigidity, the Miura origami metamaterial does not deform in shear. Consequently the in-plane shear modulus $G_{xy}$ is unbounded.

E.2 Bar-and-hinge model

The bar-and-hinge model approximates the origami folding process with a system of pin-jointed trusses and torsional springs (Schenk and Guest, 2011). By incorporating crease folding, panel bending, stretching, and shearing, the model efficiently describes essential deformation modes of an origami metamaterial (Filipov et al., 2017). As shown in Figure 17, the model depicts folding and bending with rotation along the crease and diagonal bars respectively, while stretching and shearing are equivalent to a combination of bar extension and compression. Considering axial stiffness of the bars and torsional stiffness for folding and bending, we can numerically compute the elastic properties of the metamaterial. Following the procedure introduced in (Vasudevan and Pratapa, 2024), the first step is to identify the force-displacement relation

[TABLE]

where $\mathbf{u}$ is the nodal displacement vector consisting of three translational DoF of each node, $\mathbf{f}$ collects the external applied forces, and the stiffness matrix $\mathbf{K}$ is constructed from the geometry and mechanical properties of bars and torsional elastic hinges (Filipov et al., 2017):

[TABLE]

The entries of $\mathbf{K}$ include the compatibility matrix $\mathbf{C}$ of the bar framework, the Jacobian matrices $\mathbf{J}_{B}$, $\mathbf{J}_{F}$ identifying the bending and folding hinges, and stiffness parameters of the bars and hinges $\mathbf{D}_{S}$, $\mathbf{D}_{B}$, $\mathbf{D}_{F}$, all of which detailed in (Schenk and Guest, 2011; Filipov et al., 2017). The second step is to solve Equation E.3 with appropriate loads and boundary conditions applied on the fully detailed origami metamaterial under investigation (Vasudevan and Pratapa, 2024). The bar-and-hinge model is a tessellation of Miura pattern with 20 unit cells on each edge. The loads and boundary conditions correspond to uniaxial extensions in the two in-plane directions $x$ and $y$, and an in-plane shear test. Finally, obtain the effective elastic constants from the following definitions (Vasudevan and Pratapa, 2024):

[TABLE]

where $E_{x},E_{y}$ are the in-plane Young’s modulus, $G_{xy}$ is the in-plane shear modulus, $f_{x},f_{y},f_{xy}$ are applied total loads on the boundary, $A_{x},A_{y}$ are areas of cross sections perpendicular to $x,y$ axes, and $e_{x},e_{y},e_{xy}$ are strains in their respective axial and shear directions.

E.3 Plate FE model

The plate FE method is a widely-adopted modeling strategy for thin-shelled structures including origami metamaterials. It is a preferred method in engineering design when the detailed distribution of stress and strain is required. For example, in aerospace or medical applications, engineers need such information from an FE model to minimize unpredicted failure during operation. Another advantage is that the plate FE does not involve as many assumptions on origami deformation modes as other methods. Take crease folding as an example, the rigid-panel and the bar-and-hinge models require that creases remain straight during folding, while the plate FE model allows them to curve arbitrarily. However, modeling in such details inevitably comes with significant computational costs. In many cases it becomes impractical to meet the requirements on time and resources.

Here we adopt the plate FE model introduced in (Filipov et al., 2015), with a unit cell shown in Figure 17. Panels are represented with plate FE, creases are represented with torsional springs, and the origami metamaterial is an assembly of both types of elements. The nodal DoF of a plate element at node $i$ is $\mathbf{u}_{i}=[u\;v\;w\;\theta_{x}\;\theta_{y}]^{T}$, including translational displacements $u,v$, deflection $w$, and rotations $\theta_{x},\theta_{y}$. The element stiffness matrix is

[TABLE]

where $V$ is the element volume, $\mathbf{B}$ is the matrix specifying strain-displacement relation, and $\mathbf{D}$ is the matrix of elastic coefficients for an element under plane stress. A crease is represented by a row of torsional spring elements, each connecting a pair of overlapping nodes from the two adjacent panels. The torsional stiffness of the springs sums up to the crease stiffness. The nodal DoF of a torsional spring element is $\mathbf{u}_{i}=[u\;v\;w\;\theta_{x}\;\theta_{y}\;\theta_{z}]^{T}$, including translational displacements in three directions $u,v,w$, and rotations $\theta_{x},\theta_{y},\theta_{z}$. The element stiffness matrix is

[TABLE]

where $k_{\text{lg}}$ are large values of translational and rotational stiffness compared to adjacent elements, and $K_{\textup{cr}}/n$ is the crease stiffness averaged across $n$ spring elements along the crease. In addition, to prevent detachment and tearing at the creases, we add kinematic constraints on the relative nodal displacements of each torsional spring. Each set of nodes $a,b$ of a spring element is subject to translational constraints $u_{a}-u_{b}=0$, $v_{a}-v_{b}=0$, $w_{a}-w_{b}=0$, and rotational constraints $\theta_{ya}-\theta_{yb}=0$, $\theta_{za}-\theta_{zb}=0$. The only unconstrained DoF $\theta_{x}$ lies on its local $x$ axis that aligns with the crease, meaning the torsional spring can only fold along the crease instead of detaching in any other directions. The constraints can be applied using Lagrange multipliers, for example. Finally we transform all element stiffness matrices and constraints to the global coordinates and assemble them to obtain the global stiffness matrix $\mathbf{K}$. Similar to the procedure in Section E.2, we solve Equation E.3 under three sets of loads and boundary conditions. We implement the above procedure in the commercial FE software ABAQUS, where panels are meshed with 3-node triangular plate/shell elements (S3) and connected at the creases using hinge connector elements (CONN3D2). A mesh size of $a/16$ was selected after conducting a mesh convergence test. Since we are interested in the linear elastic response, we run the simulation with a general static solver without geometric nonlinearity (Nlgeom is off). The result of displacement and reaction forces give the effective elastic constants as detailed in Section 3.

F Experimental setup

The specimen used for the compression test in Figure 12 is a 7 by 7 Miura tessellation. It is fabricated from cardboard via laser cutting (CM1290 laser cutter, SignCut Inc.) and folded by hand. Figure F.6 shows the cut pattern. Table F.1 summarizes its geometric and material properties. The geometric properties are directly measured from the specimen, while the material properties are cited from the literature (Mirzajanzadeh and Pasini, 2025).

The compression test is performed with a STEP Lab electrodynamic tester (STEP Engineering S.r.l., Resana, Treviso, Italy) using a \qty100N load cell (AEP Transducers, Congnento, Italy) and two aluminum plates as holders. The test is displacement-controlled with a strain rate of $\qty{5.26e-3}{\per}$ and a sampling frequency of \qty50Hz. We conducted 6 compression tests under the same rate.

Next, we designed a qualitative experiment to verify the homogenized results in Figure 13, where the left boundary of the metamaterial is fixed, and the right boundary has an applied load pointing [math] up to the right. A 7 by 7 Miura tessellation is laser cut from a cardboard and folded from flat. The geometric parameters of the flat pattern is identical to those in Table F.1. The left boundary of the specimen is fully fixed by gluing to a wooden panel. The right boundary is cut with an array of circular holes \qty2mm in diameter. A straight metal wire passes through the holes to act as a slider, so that the right boundary of the specimen can fold freely in-plane. A cotton string attached to slider applies an external load that evenly distribute along the slider. The load points [math] upward between $x_{1}$ and $x_{2}$ axes. Figure 13 shows the experimental setup of the fixed and slider boundary conditions, with the specimen at rest. Figure 13 shows the direction of applied load on the slider, along with the deformed specimen.

G Effective properties of curve crease origami

Curve crease origami is a typical non-rigid origami which applies to our framework. We select a pattern from the literature (Liu and James, 2024) with translational symmetry. Figure G.7 shows the origami tessellation. Figure G.7 shows the details of its unit cell. Each curved panel of the unit cell has two straight edges and two curved edges. The curved edges are defined as splines passing through a front point, a midpoint, and an end point. Table G.2 summarizes the geometry and base material properties. We report the effective compliance matrix below. Small numerical values in the matrix have been simplified to 0 because they are 2 to 4 orders of magnitude smaller than other entries, reflecting the planes of symmetry and hence anisotropy of the unit cells. The effective compliance matrix using asymptotic homogenization is

[TABLE]

The effective compliance matrix using energy-based homogenization is

[TABLE]

The effective compliance matrix of the curve crease origami metamaterial is distinct from that of Miura in Equation 47. This is because the curve crease unit cell has mirror symmetry, with one plane of symmetry parallel to the $x_{2}-x_{3}$ plane. The non-zero terms on the top right and the bottom left sections of the compliance matrix indicate coupling between in-plane strain (stress) and out-of-plane bending moments (curvature).