# Homogenization framework for rigid and non-rigid foldable origami metamaterials

Xuwen Li, Amin Jamalimehr, Mathias Legrand, Damiano Pasini

► To cite this version:

Xuwen Li, Amin Jamalimehr, Mathias Legrand, Damiano Pasini. Homogenization framework for rigid and non-rigid foldable origami metamaterials. 2025. hal-05237249

HAL Id: hal-05237249

https://hal.science/hal-05237249v1

Preprint submitted on 8 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Homogenization framework for rigid and non-rigid foldable origami metamaterials

Xuwen Li, Amin Jamalimehr, Mathias Legrand, Damiano Pasini

Department of Mechanical Engineering, McGill University, 845 Sherbrooke Street
West, Montreal, H3A0G4, Quebec, Canada

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

**Keywords:** homogenization framework, effective elastic properties, origami metamaterials, Kirchhoff-Love plate, Finite-Element modeling, rigid foldable origami, non-rigid foldable origami

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

<figure><img src="image_2.png" /><figcaption>Figure 1: Origami metamaterials. Left is a Miura origami sheet made of tape and Mylar. Right is a waterbomb origami sheet made of cardboard.</figcaption></figure>

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. 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 effective in simplifying the mechanical
analysis of large, periodic tessellations, such as cellular origami metamaterials (Cheung et al., 2014), 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. 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 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 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.

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. 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(a) shows in blue the physical domain $\Omega$ of an origami sheet. Figure 2(b) 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(c) 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(d) 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(d). 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

$$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) \quad (1a)$$

$$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), \quad (1b)$$

$$u_3(x_1, x_2, x_3) = u_{03}(x_1, x_2), \quad (1c)$$

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
<figure><img src="image_3.png" /><figcaption>Figure 2: Homogenization scheme of a single-layered, periodic origami sheet. (a) Detailed Miura sheet $\Omega$. (b) Homogeneous Kirchhoff-Love plate. (c) Mid-surface of homogenized plate $\Omega^0$. (d) Zoom-in of origami sheet showing the location of the neutral surface (superimposed orange plane) with respect to the detailed origami sheet.</figcaption></figure>

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

$$e_{11} = \partial_{x_1} u_{01}, \qquad e_{22} = \partial_{x_2} u_{02}, \qquad e_{12} = \frac{1}{2}(\partial_{x_2} u_{01} + \partial_{x_1} u_{02}), \qquad (2a)$$

$$\kappa_{11} = -\partial_{x_1}^2 u_{03}, \qquad \kappa_{22} = -\partial_{x_2}^2 u_{03}, \qquad \kappa_{12} = -\partial_{x_1} \partial_{x_2} u_{03}. \qquad (2b)$$

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

$$\begin{bmatrix} \mathbf{N} \\ \mathbf{M} \end{bmatrix} = \begin{bmatrix} \mathbf{A}^H & \mathbf{B}^H \\ \mathbf{B}^H & \mathbf{D}^H \end{bmatrix} \begin{bmatrix} \mathbf{e} \\ \boldsymbol{\kappa} \end{bmatrix}, \qquad (3)$$

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_{\alpha\beta\gamma\delta}^H$, $\mathbf{B}^H$ is the bending-extensional coupling stiffness coefficient matrix with components $B_{\alpha\beta\gamma\delta}^H$, and $\mathbf{D}^H$ is the bending stiffness coefficient matrix with components $D_{\alpha\beta\gamma\delta}^H$. 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_{\alpha\beta\gamma\delta}^H$, $B_{\alpha\beta\gamma\delta}^H$, $D_{\alpha\beta\gamma\delta}^H$ 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 a origami metamaterial with a general fold pattern. Rather, $A_{\alpha\beta\gamma\delta}^H$, $B_{\alpha\beta\gamma\delta}^H$, and $D_{\alpha\beta\gamma\delta}^H$ should all be computed independently (Vasudevan and Pratapa, 2024). The strain energy density of a homogeneous thin plate is

$$\bar{U}^H = \frac{1}{2} \begin{bmatrix} \mathbf{e} \\ \boldsymbol{\kappa} \end{bmatrix}^T \begin{bmatrix} \mathbf{A}^H & \mathbf{B}^H \\ \mathbf{B}^H & \mathbf{D}^H \end{bmatrix} \begin{bmatrix} \mathbf{e} \\ \boldsymbol{\kappa} \end{bmatrix}. \qquad (4)$$

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

<figure><img src="image_4.png" /><figcaption>Figure 3: Coordinate systems used in asymptotic homogenization. (a) A zoomed-in section of the origami sheet showing the macroscopic coordinates $\mathbf{x}$ and domain boundaries $\Gamma_h^+$, $\Gamma_h^-$, $\Gamma_h^{\text{lat}}$. (b) Domain $\Omega_h$ of the origami sheet. (c) A unit cell $Y$ scaled up by $1/\epsilon$ and the microscopic coordinates $\mathbf{y}$.</figcaption></figure>

Consider the representative origami sheet shown in Figures 3(a) and 3(b). $\Omega_h$ denotes the domain of the sheet. The upper, lower, and lateral surfaces are respectively $\Gamma_h^+$, $\Gamma_h^-$, $\Gamma_h^{\text{lat}}$. A small parameter $h$ denotes the domain height of the equivalent homogenized plate. 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 $\epsilon = \ell/n$. Since the origami sheet is single-layered, its height $h$ and unit cell size $\epsilon$ have the same order of magnitude. They are small parameters, meaning $h \to 0$ and $\epsilon \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 $\epsilon$. 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/\epsilon$, $y_2 = x_2/\epsilon$, and $y_3 = x_3/h$. The elastic stiffness of the material in the unit cell $C_{ijkl}$ is a function of $y_1$, $y_2$, and $y_3$. Before homogenization, the elastic stiffness of the full origami metamaterial $C_{ijkl}^{h\epsilon}(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

$$C_{ijkl}^{h\epsilon}(x_1, x_2, x_3) = \frac{1}{h^3} C_{ijkl} \left( \frac{x_1}{\epsilon}, \frac{x_2}{\epsilon}, \frac{x_3}{h} \right), \quad (5)$$

where the superscripts $h, \epsilon$ signifies 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^{\text{lat}}$. The elastic stiffness $C_{ijk\ell}^{h\epsilon}(x_1, x_2, x_3)$ becomes $C_{ijk\ell}^{\epsilon}(z_1, z_2, z_3)$, which is a function of the periodicity $\epsilon$ but not a function of the height $h$. The full problem of elasticity of the origami sheet is

$$ \partial_{z_1} \sigma_{i1}^{\epsilon} + \partial_{z_2} \sigma_{i2}^{\epsilon} + \frac{1}{h} \partial_{z_3} \sigma_{i3}^{\epsilon} = 0 \qquad \text{in } \Omega \qquad (6a) $$

$$ \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} \qquad \text{in } \Omega \qquad (6b) $$

$$ u_i^{\epsilon} = 0 \qquad \text{on } \Gamma^{\text{lat}} \qquad (6c) $$

$$ \sigma_{1j}^{\epsilon} n_j = t_1/h, \quad \sigma_{2j}^{\epsilon} n_j = t_2/h, \quad \sigma_{3j}^{\epsilon} n_j = t_3 \qquad \text{on } \Gamma^{\pm} \qquad (6d) $$

where $\sigma_{ij}^{\epsilon}$ is the stress tensor, $u_i^{\epsilon}$ 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 $\epsilon \to 0$. The displacement is represented by the asymptotic expansion

$$ 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 \qquad (7) $$

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 into

$$ \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 \qquad (8) $$

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 $\epsilon$. Grouping the terms with the same power of $\epsilon$ yields a set of equations that simplifies to two microscopic problems.

The first one is

$$ \partial_{y_j} [C_{ijpq} (\delta_{pk}\delta_{q\delta} + \frac{1}{2}(\partial_{y_q} w_p^{k\delta} + \partial_{y_p} w_q^{k\delta}))] = 0 \text{ on } Y, \quad \mathbf{w}^{k\delta} \text{ is } (y_1, y_2)\text{-periodic}, \qquad (9) $$

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

$$ \partial_{y_j} [C_{ijpq} (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}))] = 0 \text{ on } Y, \quad \mathbf{p}^{\gamma\delta} \text{ is } (y_1, y_2)\text{-periodic}, \qquad (10) $$

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

$$ 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}, \qquad (11) $$
<figure><img src="image_5.png" /><figcaption>Figure 4: Origami unit cell Y consisting of thin plates and their mid-surface $Y^0$.</figcaption></figure>

where $L_1$ and $L_2$ are the lengths of the unit cell in the $y_1$ and $y_2$ directions, 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, y'_2)$. 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

$$ \mathbf{e} = \begin{bmatrix} e_{11} & e_{12} \\ e_{12} & e_{22} \end{bmatrix}, \quad 2\boldsymbol{\mu} = \begin{bmatrix} 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{bmatrix}, \quad \boldsymbol{\kappa} = - \begin{bmatrix} \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{bmatrix} \quad (12) $$

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

$$ A_{\alpha\beta\gamma\delta}^{H} = \frac{1}{L_1 L_2} \int_{Y^0} C_{\zeta\eta\theta\lambda} \left( \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} \right) \left( \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} \right) dy'_1 dy'_2. \quad (13) $$

Its shorthand expression is

$$ A_{\alpha\beta\gamma\delta}^{H} = \frac{1}{L_1 L_2} \int_{Y^0} \mathbf{C}(\mathbf{e}^{0e,(\gamma\delta)} + (\boldsymbol{\mu}^{*e,(\gamma\delta)} + y_3\boldsymbol{\kappa}^{*e,(\gamma\delta)})) : (\mathbf{e}^{0e,(\alpha\beta)} + (\boldsymbol{\mu}^{*e,(\alpha\beta)} + y_3\boldsymbol{\kappa}^{*e,(\alpha\beta)})) \, dy', \quad (14) $$

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_{1111}^H, A_{2222}^H, A_{1212}^H$ and off-diagonal terms $A_{1122}^H, A_{1112}^H, A_{2212}^H$. To find the diagonal terms $A_{\alpha\beta\alpha\beta}^H$, 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)

$$ A_{\alpha\beta\alpha\beta}^{H} = \frac{2}{L_1 L_2} U. \quad (15) $$
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

$$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, \qquad (16)$$

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

$$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}. \qquad (17)$$

Similarly,

$$D_{\alpha\beta\gamma\delta}^H = \frac{1}{L_1 L_2} \int_Y C_{ijkl} (y_3 \delta_{k\gamma} \delta_{\ell\delta} + \frac{1}{2} (\partial_{y_l} p_k^{\gamma\delta} + \partial_{y_k} p_l^{\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})) \, dy \quad (18)$$

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

$$D_{\alpha\beta\gamma\delta}^H = \frac{1}{L_1 L_2} \int_{Y^0} \mathbf{C}(\mathbf{e}^{0\kappa, (\gamma\delta)} + (\boldsymbol{\mu}^{*\kappa, (\gamma\delta)} + y_3 \boldsymbol{\kappa}^{*\kappa, (\gamma\delta)})) : (\mathbf{e}^{0\kappa, (\alpha\beta)} + (\boldsymbol{\mu}^{*\kappa, (\alpha\beta)} + y_3 \boldsymbol{\kappa}^{*\kappa, (\alpha\beta)})) \, dy, \qquad (19)$$

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 term as

$$D_{\alpha\beta\alpha\beta}^H = \frac{2}{L_1 L_2} U. \qquad (20)$$

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

$$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}. \qquad (21)$$

Another set of effective stiffness coefficients

$$B_{\alpha\beta\gamma\delta}^H = \frac{1}{L_1 L_2} \int_Y C_{ijkl} (y_3 \delta_{k\gamma} \delta_{\ell\delta} + \frac{1}{2} (\partial_{y_l} p_k^{\gamma\delta} + \partial_{y_k} p_l^{\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})) \, dy, \quad (22)$$

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

$$B_{\alpha\beta\gamma\delta}^H = \frac{1}{L_1 L_2} \int_{Y^0} \mathbf{C}(\mathbf{e}^{0\kappa, (\gamma\delta)} + (\boldsymbol{\mu}^{*\kappa, (\gamma\delta)} + y_3 \boldsymbol{\kappa}^{*\kappa, (\gamma\delta)})) : (\mathbf{e}^{0e, (\alpha\beta)} + (\boldsymbol{\mu}^{*e, (\alpha\beta)} + y_3 \boldsymbol{\kappa}^{*e, (\alpha\beta)})) \, dy. \quad (23)$$
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

$$B_{\alpha\beta\gamma\delta}^H = \frac{1}{L_1L_2}U - \frac{A_{\alpha\beta\alpha\beta}^H}{2} - \frac{D_{\gamma\delta\gamma\delta}^H}{2}. \qquad (24)$$

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

$$\left\{
\begin{array}{ll}
\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, \quad \bar{u}_{\alpha} = 0 & \text{in } \Gamma^0,
\end{array}
\right.
\qquad (25)$$

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 ds$, $Q_\alpha$ are averaged bending moments defined as $Q_\alpha = \frac{h}{|Y|} \int_{\Gamma^\pm} y_3 t_\alpha^\pm ds$, $\bar{u}_\alpha$ are the homogenized displacements, and $u_3^{(0)}$ 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 at https://github.com/Wendyli119/OMhomogenization.git. 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(c). This choice simplifies the specification of boundary conditions in the microscopic problem.
<table><caption>Table 1: Geometric and material parameters of Miura origami metamaterial.</caption><tr><td>Panel side length</td><td>a (mm)</td><td>20</td></tr><tr><td>Panel side length</td><td>b (mm)</td><td>20</td></tr><tr><td>Sector angle</td><td>γ (°)</td><td>60</td></tr><tr><td>Initial fold angle</td><td>θ<sub>0</sub> (°)</td><td>30</td></tr><tr><td>Panel thickness</td><td>t (mm)</td><td>0.13</td></tr><tr><td>Unit length crease stiffness</td><td>K<sub>cr</sub>/a (N)</td><td>0.1</td></tr><tr><td>Base material Young's modulus</td><td>E (MPa)</td><td>4000</td></tr><tr><td>Base material Poisson's ratio</td><td>ν</td><td>0.38</td></tr></table>

<figure><img src="image_6.png" /><figcaption>Figure 5: Miura unit cells and geometric parameters. (a) A Miura tessellation. (b) A common choice of Miura unit cell. (c) Equivalent unit cell adopted in this work.</figcaption></figure>

Discretizing the microscopic problems exposed in Equations (9) and (10) leads to solving two systems of linear equations $K\mathbf{w}^{\gamma\delta} = \mathbf{f}^0$ and $K\mathbf{q}^{\gamma\delta} = \mathbf{f}^0$, where $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 0 everywhere in the unit cell. $\mathbf{f}^0$ is expressed as $\mathbf{f}^0 = \int_Y \mathbf{B}^T\mathbf{E}\boldsymbol{\chi}^0 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

$$ \begin{align} e_{11} = 1, \quad \boldsymbol{\chi}_{\text{node}}^0 &= \begin{bmatrix} y_1 \\ 0 \\ 0 \end{bmatrix}; \quad e_{22} = 1, \quad \boldsymbol{\chi}_{\text{node}}^0 = \begin{bmatrix} 0 \\ y_2 \\ 0 \end{bmatrix}; \quad e_{12} = 1, \quad \boldsymbol{\chi}_{\text{node}}^0 = \frac{1}{2} \begin{bmatrix} y_2 \\ y_1 \\ 0 \end{bmatrix} \nonumber \\ \kappa_{11} = 1, \quad \boldsymbol{\chi}_{\text{node}}^0 &= \frac{1}{2} \begin{bmatrix} 2y_1y_3 \\ 0 \\ -y_1^2 \end{bmatrix}; \quad \kappa_{22} = 1, \quad \boldsymbol{\chi}_{\text{node}}^0 = \frac{1}{2} \begin{bmatrix} 0 \\ 2y_2y_3 \\ -y_2^2 \end{bmatrix}; \quad \kappa_{12} = 1, \quad \boldsymbol{\chi}_{\text{node}}^0 = \frac{1}{2} \begin{bmatrix} y_2y_3 \\ y_1y_3 \\ -y_1y_2 \end{bmatrix} \tag{26} \end{align} $$

where $\boldsymbol{\chi}_{\text{node}}^0$ 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.

<figure><img src="image_7.png" /><figcaption>Figure 6: Step 1 of asymptotic homogenization showing six of the deformation modes with unit strain or unit curvature. To ensure uniform strain or curvature distribution, displacements are prescribed nodally as shown in the inset. Displacements are scaled.</figcaption></figure>

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

$$e_{11} = 1, \quad e_{22} = 1, \quad \chi_{\text{node}}^{0} = \begin{bmatrix} y_1 \\ 0 \\ 0 \end{bmatrix} + \begin{bmatrix} 0 \\ y_2 \\ 0 \end{bmatrix} = \begin{bmatrix} y_1 \\ y_2 \\ 0 \end{bmatrix}. \qquad (27)$$

The initial strain loading $f^0$ is then found from

$$K \chi^0 = f^0. \qquad (28)$$

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 $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

$$K \chi^* = f^0, \qquad (29)$$

where $f^0$ is the initial strain loading obtained from the previous step, and $\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 the $(y_1, y_2)$ directions, 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

$$ \chi^{*(l)} - \chi^{*(r)} = 0, \qquad (30a) $$

$$ \chi^{*(t)} - \chi^{*(b)} = 0, \qquad (30b) $$

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

The initial strain loadings $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 $\chi^*$. Figure 7 shows $\chi^*$ corresponding to the deformation mode $e_{11} = 1$.

<figure><img src="image_8.png" /><figcaption>Figure 7: Step 2 of asymptotic homogenization solving the microscopic problem. (a) Periodic boundary conditions are applied on the boundary displacements according to Equation (30). (b) Numerical solutions of the characteristic displacement field $\chi^*$. Displacements are scaled.</figcaption></figure>

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

$$ \tilde{\chi} = \chi^0 - \chi^*. \qquad (31) $$

$\tilde{\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{\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).
<figure><img src="image_9.png" /><figcaption>Figure 8: Step 3 of asymptotic homogenization computing unit cell strain energy $U$. (a) Difference between the unit strain displacement field $\chi^0$ and the characteristic displacement field $\chi^*$. Only 6 of the deformation modes are shown with lighting effects in 3D rendering. Displacements are scaled. (b) Unit cell strain energy $U$. Colors represent the strain energy magnitudes of the corresponding deformation modes.</figcaption></figure>

The final step is the homogenized problem. The origami sheet is simplified into a thin plate. The homogenized plate $\Omega$, shown in Figure 2(c), 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 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 material 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)

$$u_{\alpha}^{(r)} = u_{\alpha}^{(l)} + e_{\alpha\beta}L_{\beta}, \qquad (32a)$$

$$u_{3}^{(r)} = u_{3}^{(l)}, \qquad (32b)$$

$$u_{\alpha}^{(t)} = u_{\alpha}^{(b)} + e_{\alpha\beta}L_{\beta}, \qquad (32c)$$

$$u_{3}^{(t)} = u_{3}^{(b)}, \qquad (32d)$$

where $u_{\alpha}^{(r)}$, $u_{\alpha}^{(l)}$, $u_{\alpha}^{(t)}$, and $u_{\alpha}^{(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 in $y_{\beta}$ directions. 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)

$$\begin{bmatrix} u_{1}^{(r)} - u_{1}^{(l)} \\ u_{2}^{(r)} - u_{2}^{(l)} \\ u_{3}^{(r)} - u_{3}^{(l)} \end{bmatrix} = \frac{1}{2} \begin{bmatrix} (x_{1}^{(r)} - x_{1}^{(l)})(h - 2x_{3}^{(r)}) & 0 & 0 \\ 0 & 0 & (x_{1}^{(r)} - x_{1}^{(l)})(h - 2x_{3}^{(r)})/2 \\ (x_{1}^{(r)} - x_{1}^{(l)})(x_{1}^{(r)} + x_{1}^{(l)}) & 0 & x_{1}^{(r)}x_{2}^{(r)} - x_{1}^{(l)}x_{2}^{(l)} \end{bmatrix} \begin{bmatrix} \kappa_{11} \\ \kappa_{22} \\ 2\kappa_{12} \end{bmatrix}, \quad (33a)$$

$$\begin{bmatrix} u_{1}^{(t)} - u_{1}^{(b)} \\ u_{2}^{(t)} - u_{2}^{(b)} \\ u_{3}^{(t)} - u_{3}^{(b)} \end{bmatrix} = \frac{1}{2} \begin{bmatrix} 0 & 0 & (x_{2}^{(t)} - x_{2}^{(b)})(h - 2x_{3}^{(t)})/2 \\ 0 & (x_{2}^{(t)} - x_{2}^{(b)})(h - 2x_{3}^{(t)}) & 0 \\ 0 & (x_{2}^{(t)} - x_{2}^{(b)})(x_{2}^{(t)} + x_{2}^{(b)}) & x_{1}^{(t)}x_{2}^{(t)} - x_{1}^{(b)}x_{2}^{(b)} \end{bmatrix} \begin{bmatrix} \kappa_{11} \\ \kappa_{22} \\ 2\kappa_{12} \end{bmatrix}. \quad (33b)$$

where $x_{i}^{(r)}$, $x_{i}^{(l)}$, $x_{i}^{(t)}$, $x_{i}^{(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

$$\bar{U} = \frac{1}{2L_1 L_2} \int_Y \hat{\sigma}_{ij} \hat{e}_{ij} \, dy_1 \, dy_2 \qquad (34)$$

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

$$\bar{U}^H = \bar{U}, \qquad (35)$$

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

$$u_1^{(r)} - u_1^{(l)} - (y_1^{(r)} - y_1^{(l)}) = 0, \qquad (36a)$$

$$u_2^{(r)} - u_2^{(l)} = 0, \qquad (36b)$$

$$u_3^{(r)} - u_3^{(l)} = 0, \qquad (36c)$$

where $u_i^{(r)}$, $u_i^{(l)}$ are the right and left nodal displacements in $y_i$, and $y_1^{(r)}$, $y_1^{(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

$$u_i^{(t)} - u_i^{(b)} = 0. \qquad (37)$$

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

$$e_{11} = 1, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \begin{bmatrix} y_1^{(r)} - y_1^{(l)} \\ 0 \\ 0 \end{bmatrix}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \mathbf{0}, \qquad (38)$$

$$e_{22} = 1, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \mathbf{0}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \begin{bmatrix} 0 \\ y_2^{(t)} - y_2^{(b)} \\ 0 \end{bmatrix} \qquad (39)$$

$$e_{12} = 2, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \begin{bmatrix} 0 \\ y_1^{(t)} - y_1^{(b)} \\ 0 \end{bmatrix}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \begin{bmatrix} y_2^{(r)} - y_2^{(l)} \\ 0 \\ 0 \end{bmatrix}. \qquad (40)$$

For out-of-plane deformation modes,

$$\kappa_{11} = 1, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \begin{bmatrix} \left(\frac{h}{2} - y_3^{(r)}\right)(y_1^{(r)} - y_1^{(l)}) \\ 0 \\ \frac{1}{2}(y_1^{(r)}{}^2 - y_1^{(l)}{}^2) \end{bmatrix}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \mathbf{0}, \qquad (41)$$

$$\kappa_{22} = 1, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \mathbf{0}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \begin{bmatrix} 0 \\ \left(\frac{h}{2} - y_3^{(t)}\right)(y_2^{(t)} - y_2^{(b)}) \\ \frac{1}{2}(y_2^{(t)}{}^2 - y_2^{(b)}{}^2) \end{bmatrix}, \qquad (42)$$

$$\kappa_{12} = 2, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \begin{bmatrix} 0 \\ \frac{1}{2}(\frac{h}{2} - y_3^{(r)})(y_1^{(r)} - y_1^{(l)}) \\ \frac{1}{2}(y_1^{(r)} y_2^{(r)} - y_1^{(l)} y_2^{(l)}) \end{bmatrix}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \begin{bmatrix} \frac{1}{2}(\frac{h}{2} - y_3^{(t)})(y_2^{(t)} - y_2^{(b)}) \\ 0 \\ \frac{1}{2}(y_1^{(t)} y_2^{(t)} - y_1^{(b)} y_2^{(b)}) \end{bmatrix} \qquad (43)$$

where $y_3^{(r)} = y_3^{(l)}$ and $y_3^{(t)} = y_3^{(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.

<figure><img src="image_10.png" /><figcaption>Figure 9: Step 1 of energy-based homogenization identifying deformation modes. (a) Four boundaries of the unit cell. (b) Deformed boundaries after applying periodic boundary conditions. The overall strain is $e_{11} = 1$. (c) Six of the deformation modes under periodic boundary conditions with lighting effects in 3D rendering. Displacements are scaled.</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>Figure 10: Step 2 of energy-based homogenization computes the total unit cell strain energy U of each deformation mode in step 1. Colors represent the strain energy magnitudes of the corresponding deformation modes.</figcaption></figure>

In addition, we write the 15 linear combinations of the periodic boundary conditions in Equations (38) to (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

$$e_{11} = e_{22} = 1, \quad \mathbf{u}^{(r)} - \mathbf{u}^{(l)} = \begin{bmatrix} y_1^{(r)} - y_1^{(l)} \\ 0 \\ 0 \end{bmatrix}, \quad \mathbf{u}^{(t)} - \mathbf{u}^{(b)} = \begin{bmatrix} 0 \\ y_2^{(t)} - y_2^{(b)} \\ 0 \end{bmatrix}. \quad (44)$$

Similar to Equations (38) to (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_{ijkl}^H$, $B_{ijkl}^H$, and $D_{ijkl}^H$ from Equation (35) (Wang et al., 2009). Each strain energy $U$ computed from Equations (38) to (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) to (43), corresponds to a total strain energy in the unit cell $U$. From the energy equivalence, we have

$$A_{\alpha\beta\gamma\delta}^{H} = \frac{2U}{L_{1}L_{2}e_{\alpha\beta}e_{\gamma\delta}} \quad \text{and} \quad D_{\alpha\beta\gamma\delta}^{H} = \frac{2U}{L_{1}L_{2}\kappa_{\alpha\beta}\kappa_{\gamma\delta}}, \quad (45)$$

where $L_1$ and $L_2$ are the lengths of the unit cell in $y_1$ and $y_2$ directions, 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_{1111}^H$ describes the coupling effect between uniaxial extension in $e_{11}$ and bending $\kappa_{11}$. Energy equivalence indicates that $B_{1111}^H$ 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

$$\begin{align}
A_{\alpha\beta\gamma\delta}^{H} &= \frac{1}{e_{\alpha\beta}e_{\gamma\delta}} \left( \frac{U}{L_1 L_2} - \frac{1}{2e_{\alpha\beta}e_{\alpha\beta}} A_{\alpha\beta\alpha\beta}^{H} - \frac{1}{2e_{\gamma\delta}e_{\gamma\delta}} A_{\gamma\delta\gamma\delta}^{H} \right), \nonumber \\
B_{\alpha\beta\gamma\delta}^{H} &= \frac{1}{e_{\alpha\beta}\kappa_{\gamma\delta}} \left( \frac{U}{L_1 L_2} - \frac{1}{2e_{\alpha\beta}e_{\alpha\beta}} A_{\alpha\beta\alpha\beta}^{H} - \frac{1}{2\kappa_{\gamma\delta}\kappa_{\gamma\delta}} D_{\gamma\delta\gamma\delta}^{H} \right), \tag{46} \\
D_{\alpha\beta\gamma\delta}^{H} &= \frac{1}{\kappa_{\alpha\beta}\kappa_{\gamma\delta}} \left( \frac{U}{L_1 L_2} - \frac{1}{2\kappa_{\alpha\beta}\kappa_{\alpha\beta}} D_{\alpha\beta\alpha\beta}^{H} - \frac{1}{2\kappa_{\gamma\delta}\kappa_{\gamma\delta}} D_{\gamma\delta\gamma\delta}^{H} \right). \nonumber
\end{align}$$

Note that the results from uniaxial extension, shear, bending, and twisting are required for $A_{\alpha\beta\alpha\beta}^H$ and $D_{\alpha\beta\alpha\beta}^H$ 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° 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):

$$ \begin{bmatrix} e_{11} \\ e_{22} \\ 2e_{12} \\ \kappa_{11} \\ \kappa_{22} \\ 2\kappa_{12} \end{bmatrix} = \begin{bmatrix} \frac{1}{E_1} & -\frac{\nu_{21}^s}{E_2} & 0 & 0 & 0 & 0 \\ -\frac{\nu_{12}^s}{E_1} & \frac{1}{E_2} & 0 & 0 & 0 & 0 \\ 0 & 0 & \frac{1}{G_{12}} & 0 & 0 & 0 \\ 0 & 0 & 0 & \frac{1}{M_1} & -\frac{\nu_{21}^b}{M_2} & 0 \\ 0 & 0 & 0 & -\frac{\nu_{12}^b}{M_1} & \frac{1}{M_2} & 0 \\ 0 & 0 & 0 & 0 & 0 & \frac{1}{T_{12}} \end{bmatrix} \begin{bmatrix} \sigma_{11} \\ \sigma_{22} \\ \sigma_{12} \\ \mu_{11} \\ \mu_{22} \\ \mu_{12} \end{bmatrix}, \quad (47) $$

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_{12}^b$, bending moduli $M_1$ and $M_2$, and twisting modulus $T_{12}$. These relations, namely the effective elastic constants, can be determined from the effective stiffness coefficients as

$$ \mathbf{S} = h[\mathbf{ABD}]^{-1} \quad (48) $$

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

Figure 11(a) 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(b) 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(c) 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(d) 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
<figure><img src="image_12.png" /><figcaption>Figure 11: Detailed plate FE model under mechanical tests: (a,b) uniaxial extension, (c) in-plane shear, (d,e) bending, and (f) twisting.</figcaption></figure>

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_{12}^b = -\kappa_2/\kappa_1$. Figure 11(e) 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(f) 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.

Appendix 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 3.47%, the energy-based plate FE 3.20%, the bar-and-hinge 3.83%, and the rigid-panel 32.1%. The results of $G_{12}$ and $nu_{12}^b$ 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 11.1%, energy-based plate FE 12.9%,) than the bar-and-hinge model (error 161%).
## 4. Results

We can now apply both homogenization methods to determine the effective in-plane and out-of-plane properties of Miura origami sheets with geometric and material parameters specified in Table 1. Appendix C provides the results in matrix form. 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 12(a) shows a homogenized plate Ω with equal dimensions as the tessellation. The above compliance matrices are applied to the homogenized models as detailed in the last steps of Section 2.2.2 and Section 2.3.2. There is a fixed boundary condition on the left boundary and a load of 25 N evenly distributed on the right boundary, pointing 45° to the upper right. Figures 12(b) and 12(c) presents the resulting displacement fields in the homogenized models. Identical load and boundary condition are applied to the detailed plate FE model as the baseline, whose displacement is shown in Figure 12(d). 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.

<figure><img src="image_13.png" /><figcaption>Figure 12: Comparison between homogenized and detailed displacement magnitudes. (a) Boundary condition and load chosen for demonstration. (b) Deformed asymptotic homogenization model. (c) Deformed energy-based homogenization model. Both homogenized models (b) and (c) have 575 elements. (d) Detailed model with 317715 elements.</figcaption></figure>

### 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 $θ_0$, as shown in Figure 13. It 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 $θ_0$ as the varying parameter due to its strong influence on the folding kinematics. It also allows direct comparison with results in the literature. In this parametric study, $θ_0$ varies from 1° to 80°. 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_{12}^b$, bending moduli $M_1, M_2$, and twisting modulus $T_{12}$.

<figure><img src="image_14.png" /><figcaption>Figure 13: Two configurations of Miura unit cell with distinct fold angles $\theta_0$.</figcaption></figure>

Figures 14(a) and 14(e) 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$, 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 14(b) and 14(c). Figure 14(d) 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 14(f) 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 14(g) and 14(h), which is qualitatively in agreement with the response of a truss-based lattice (Vasudevan and Pratapa, 2024).
<figure><img src="image_15.png" /><figcaption>Figure 14: In-plane and out-of-plane effective elastic constants of Miura origami metamaterial. Effective Young's moduli $E_1$, $E_2$, shear modulus $G_{12}$, bending moduli $M_1$, $M_2$, and twisting modulus $T_{12}$ are normalized against the base material Young's modulus $E_{\text{base}}$. Insets show corresponding deformation modes of a detailed FE model. A subset of data points is shown for visual clarity. (★) Asymptotic plate FE; (▲) Energy-based plate FE; (●) Detailed plate FE.</figcaption></figure>
## 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 15 summarizes all modeling methods under comparison. As with the plate FE results presented in Section 4.1, the

<figure><img src="image_16.png" /><figcaption>Figure 15: Comparison of modeling methods applied to a Miura unit cell for demonstrative purposes.</figcaption></figure>

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_{12}^b$, $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 16(a), 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, 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_{12}^b$ in Figure 16(e) is close to the negative of $\nu_{12}$ in Figure 16(a).
<figure><img src="image_17.png" /><figcaption>Figure 16: Comparison of the homogenization framework proposed in this work with methods in the literature. Sub-plots show the variation of effective elastic constants over the fold angles $θ_0$ between 0° and 80°. Except for the Poisson's ratios $v_{12}$ and $v_{12}^b$, the effective elastic constants are normalized against the base material Young's modulus $E_{\text{base}}$. Insets show zoomed-in regions of the plots. A subset of data points is shown for visual clarity. (★) Asymptotic plate FE; (▲) Energy-based plate FE; (●) Detailed plate FE; (♦) Homogenized/detailed bar-and-hinge; (---) Rigid panel.</figcaption></figure>
Figures 16(b) and 16(c) 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 16(f), 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 16(g), $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 16(h). 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 1.3 N-mm for each 20 mm long crease. To put numbers in context, the crease stiffness $K_{cr}$ of a typical origami metamaterial is smaller than or similar to the panel bending stiffness $K_p$. The panel bending stiffness is expressed as

$$K_p = \left(0.55 - 0.42 \frac{\sum \alpha}{\pi}\right) \frac{E_{\text{base}} t^3}{12(1 - \nu^2)} \left(\frac{D_s}{t}\right)^{1/3}, \quad (49)$$

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 17 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_{cr}/K_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 17(a) and 17(e) shows that the crease stiffness has negligible influence on both in-plane and bending Poisson's ratios $\nu_{12}, \nu_{12}^b$. This suggests that when the crease stiffness is smaller than or equal to the panel bending stiffness,
<figure><img src="image_18.png" /><figcaption>Figure 17: Variation of effective elastic constants with respect to crease stiffness. Effective Young's moduli $E_1$, $E_2$, shear modulus $G_{12}$, bending moduli $M_1$, $M_2$, and twisting modulus $T_{12}$ are normalized against the base material Young's modulus $E_{\text{base}}$. Crease stiffness $K_{\text{cr}}$ is normalized against panel bending stiffness $K_p$. A subset of data points is shown for visual clarity. (★) Asymptotic plate FE; (▲) Energy-based plate FE; (●) Detailed plate FE; (■) Homogenized/detailed bar-and-hinge; (---) Rigid panel.</figcaption></figure>
the folding kinematics of Miura origami is consistent. In Figures 17(b) and 17(c), 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 17(d) 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 17(f) and 17(g), meaning the stiffer the creases, the more difficult the origami sheet bends. The twisting modulus in Figure 17(h) 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 $v_{12}$, $v_{12}^b$, $G_{12}$, $M_1$, and $M_2$ show excellent agreement. Figure 17(d) 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 17(b) and 17(c). 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_{cr}/K_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 17(h) 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 16 and 17 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 18(a) 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 deformation modes such as the emergence of curvature along straight creases, local strain at vertices and warpage in panels.

<figure><img src="image_19.png" /><figcaption>Figure 18: Unit cell deformation modes in various simulation methods. (a) Miura unit cell with rigid panels has 1 DoF. (b) Stretched origami sheet under lateral constraints. Deformed unit cells at the center of the stretched sheet modeled with (c) plate FE and (d) bar-and-hinge plotted with identical scale factor. (e) Twisted origami sheet. Deformed unit cells at the center of the twisted sheet modeled with (f) plate FE and (g) bar-and-hinge plotted with identical scale factor.</figcaption></figure>

## 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 18(b) 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 18(c) 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 18(d) 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 18(e) 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 18(f) and 18(g). 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 18(g) 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 compare 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 8% 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_{\alpha\beta\gamma\delta}^H$ 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 (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: https://github.com/Wendyli119/OMhomogenization.
git

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 discus-
sions 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 0.05 N. We increase the number of elements along each
side of the rhombic panel and plot the variation of x₁-direction Young's modulus E₁. 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°, 80°), all models share identical geometry
<figure><img src="image_20.png" /><figcaption>Figure A.1: Mesh convergence test and validation with detailed FE model. (a) In-plane and (b) out-of-plane effective properties from energy-based homogenization. (c) In-plane and (d) out-of-plane effective properties from asymptotic homogenization.</figcaption></figure>

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

$$ \text{Error}(\theta_0) = \frac{|P(\text{model}) - P(\text{detailed plate FE})|}{|P(\text{detailed plate FE})|} \times 100\%, \quad (B.1) $$

where P represents any of the elastic constants $\nu_{12}$, $E_1$, $E_2$, $\nu_{12}^b$, $M_1$, $M_2$, and $T_{12}$. Their errors under varying initial fold angles $\theta_0$ are presented in Figure B.3. The error over the fold angles is averaged as

$$ \text{Averaged error} = \sqrt{\frac{1}{\theta_0^{\text{max}} - \theta_0^{\text{min}}} \int_{\theta_0^{\text{min}}}^{\theta_0^{\text{max}}} \text{Error}^2(\theta_0) \, d\theta_0}, \quad (B.2) $$

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

$$ \text{Maximum error} = \max_{\theta_0 \in (0^{\circ}, 80^{\circ})} \text{Error}(\theta_0). \quad (B.3) $$

Errors of all elastic constants are compared in Figure B.2. On average, Poisson’s ratios $\nu_{12}$,
<figure><img src="image_21.png" /><figcaption>Figure B.2: Comparison of errors between the proposed homogenization framework and the literature. (a) Averaged error and (b) maximum error of in-plane and out-of-plane elastic constants. Both plots describe errors within a range of fold angles between 0° and 80°, while other parameters are identical to those in Table 1. (■) Asymptotic plate FE; (■) Energy-based plate FE; (■) Homogenized/detailed bar-and-hinge; (■) Rigid panel.</figcaption></figure>

$\nu_{12}^b$ 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 53%.

### 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). The compliance matrix from asymptotic homogenization reads

$$ \begin{bmatrix} 10.6 & 23.8 & 0.000245 & 0.0589 & -0.133 & 0.0369 \\ 23.8 & 54.2 & 0.000544 & 0.133 & -0.300 & 0.0847 \\ 0.000245 & 0.000544 & 0.0729 & 0.00000861 & -0.0000201 & -0.000105 \\ 0.0589 & 0.133 & 0.00000861 & 0.787 & -1.76 & -0.0377 \\ -0.133 & -0.300 & -0.0000201 & -1.76 & 4.02 & 0.0865 \\ 0.0369 & 0.0847 & -0.000105 & -0.0377 & 0.0865 & 22.8 \end{bmatrix} \quad (C.1) $$

while for energy-based homogenization,

$$ \begin{bmatrix} 10.8 & 24.2 & 0.000444 & -0.00493 & 0.0102 & -0.179 \\ 24.2 & 55.1 & 0.00114 & -0.0103 & 0.0211 & -0.405 \\ 0.000444 & 0.00114 & 0.0732 & -0.0000178 & 0.00000321 & -0.000449 \\ -0.00493 & -0.0103 & -0.0000178 & 0.817 & -1.83 & -0.0726 \\ 0.0102 & 0.0211 & 0.00000321 & -1.83 & 4.17 & 0.165 \\ -0.179 & -0.405 & -0.000449 & -0.0726 & 0.165 & 23.98 \end{bmatrix}. \quad (C.2) $$
<figure><img src="image_22.png" /><figcaption>Figure B.3: Errors of homogenization and theoretical models compared to the detailed plate FE model under various initial fold angles $\theta_0$. (★) Asymptotic plate FE; (▲) Energy-based plate FE; (■) Homogenized/detailed bar-and-hinge; (---) Rigid panel.</figcaption></figure>

The models of both matrices use geometric and material parameters in Table 1. The units in the matrix are MPa⁻¹ for in-plane terms, mm/N for in-plane and out-of-plane coupling terms, and N⁻¹ for out-of-plane terms.
## D. Flowcharts of homogenization framework

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(a) 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_{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. Appendices E.1 to 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):

$$ \nu_{xy} = -\frac{e_y}{e_x} = -(\cos \theta \tan \gamma)^2 \qquad (E.1) $$

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

$$ E_x = \frac{4K_{cr}((1 - (\sin \gamma \sin(\rho/2)^2) + \cos^2 \gamma))}{(1 - (\sin \gamma \sin(\rho/2)^2) \sin \gamma \cos(\rho/2))}, \qquad (E.2a) $$

$$ E_y = \frac{4K_{cr}((1 - (\sin \gamma \sin(\rho/2)^2) + \cos^2 \gamma))}{(1 - (\sin \gamma \sin(\rho/2)^2)^{1/2} \cos \gamma \sin^2 \gamma \sin \rho)}, \qquad (E.2b) $$

where $\rho$, chosen as the DoF in (Wei et al., 2013), is a dihedral angle shown in Figure 2(a). Due to panel rigidity, the Miura origami metamaterial does not deform in shear. Consequently the in-plane shear modulus $G_{xy}$ is unbounded.
<figure><img src="image_23.png" /></figure>

① Displacement field equivalent to unit strain or curvature

Prescribed nodal displacements $\chi_{\text{node}}^0$

<figure><img src="image_24.png" /></figure>

$e_{11} = 1$

<figure><img src="image_25.png" /></figure>

$e_{12} = 1$

<figure><img src="image_26.png" /></figure>

$\kappa_{11} = 1$

<figure><img src="image_27.png" /></figure>

$\kappa_{12} = 1$

Given: 21 deformation modes $\chi_{\text{node}}^0$

Find: initial strain loadings $\mathbf{f}^0$

② Microscopic problem

Given:
periodic boundary conditions,
boundary conditions
preventing rigid body motion,
and initial strain loading $\mathbf{f}^0$

<figure><img src="image_28.png" /></figure>

Find: characteristic displacement $\mathbf{\chi}^*$

<figure><img src="image_29.png" /></figure>

$\mathbf{\chi}^*, e_{11} = 1$

③ Given: difference between nodal displacements $\chi^0$ and $\chi^*$

Find: unit cell strain energy $U$

<figure><img src="image_30.png" /></figure>

④ Effective elastic constants

$A_{\alpha\beta\gamma\delta}^H, B_{\alpha\beta\gamma\delta}^H, D_{\alpha\beta\gamma\delta}^H$

<figure><img src="image_31.png" /><figcaption>Figure D.4: Asymptotic homogenization flowchart</figcaption></figure>
<figure><img src="image_32.png" /></figure>

① Given: 21 periodic boundary conditions

<figure><img src="image_33.png" /></figure>

Example of first 6 deformation modes

<figure><img src="image_34.png" /></figure>

Find: unit cell strain energy U

<figure><img src="image_35.png" /></figure>

② Effective elastic constants

<figure><img src="image_36.png" /><figcaption>Figure D.5: Energy-based homogenization flowchart</figcaption></figure>
## 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 16, 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

$$ \mathbf{K}\mathbf{u} = \mathbf{f}, \qquad (E.3) $$

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

$$ \mathbf{K} = \begin{bmatrix} \mathbf{C} \\ \mathbf{J}_B \\ \mathbf{J}_F \end{bmatrix}^T \begin{bmatrix} \mathbf{D}_S & \mathbf{0} & \mathbf{0} \\ \mathbf{0} & \mathbf{D}_B & \mathbf{0} \\ \mathbf{0} & \mathbf{0} & \mathbf{D}_F \end{bmatrix} \begin{bmatrix} \mathbf{C} \\ \mathbf{J}_B \\ \mathbf{J}_F \end{bmatrix}. \qquad (E.4) $$

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

$$ \nu_{xy} = -\frac{e_y}{e_x}, \quad E_x = \frac{f_x}{A_x e_x}, \quad E_y = \frac{f_y}{A_y e_y}, \quad G_{xy} = \frac{f_{xy}}{2A_x e_{xy}}, \qquad (E.5) $$

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

$$ \mathbf{K}^e = \iiint_V \mathbf{B}^T \mathbf{D} \mathbf{B} \, dv, \qquad (E.6) $$

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

$$ \mathbf{K}^e = \begin{bmatrix} \mathbf{K}_n & -\mathbf{K}_n \\ -\mathbf{K}_n & \mathbf{K}_n \end{bmatrix} \quad \text{with} \quad \mathbf{K}_n = \begin{bmatrix} k_{lg} & 0 & 0 & 0 & 0 & 0 \\ 0 & k_{lg} & 0 & 0 & 0 & 0 \\ 0 & 0 & k_{lg} & 0 & 0 & 0 \\ 0 & 0 & 0 & K_{cr}/n & 0 & 0 \\ 0 & 0 & 0 & 0 & k_{lg} & 0 \\ 0 & 0 & 0 & 0 & 0 & k_{lg} \end{bmatrix}, \qquad (E.7) $$

where $k_{lg}$ are large values of translational and rotational stiffness compared to adjacent elements, and $K_{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 Appendix 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.
References

Almessabi, A., Li, X., Jamalimehr, A., Pasini, D., 2024. Reprogramming multi-stable snapping and energy dissipation in origami metamaterials through panel confinement. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 382, 20240005. doi:10.1098/rsta.2024.0005.

Cai, Y., Xu, L., Cheng, G., 2014. Novel numerical implementation of asymptotic homogenization method for periodic plate structures. International Journal of Solids and Structures 51, 284-292. doi:10.1016/j.ijsolstr.2013.10.003.

Caillerie, D., Nedelec, J.C., 1984. Thin elastic and periodic plates. Mathematical Methods in the Applied Sciences 6, 159-191. doi:10.1002/mma.1670060112.

Castle, T., Sussman, D.M., Tanis, M., Kamien, R.D., 2016. Additive lattice kirigami. Science advances 2, e1601258. doi:10.1126/sciadv.1601258.

Cheng, G.D., Cai, Y.W., Xu, L., 2013. Novel implementation of homogenization method to predict effective properties of periodic materials. Acta Mechanica Sinica 29, 550-556. doi:10.1007/s10409-013-0043-0.

Cheung, K.C., Tachi, T., Calisch, S., Miura, K., 2014. Origami interleaved tube cellular materials. Smart Materials and Structures 23, 094012. doi:10.1088/0964-1726/23/9/094012.

Chi, Y., Li, Y., Zhao, Y., Hong, Y., Tang, Y., Yin, J., 2022. Bistable and multistable actuators for soft robots: Structures, materials, and functionalities. Advanced Materials 34, 2110384. doi:10.1002/adma.202110384.

Dang, X., Chen, S., Acha, A.E., Wu, L., Pasini, D., 2025. Shape and topology morphing of closed surfaces integrating origami and kirigami. Science Advances 11, eads5659. doi:10.1126/sciadv.ads5659.

Eskandari, S., Shahryari, B., Akbarzadeh, A., 2024. Unravelling size-dependent and coupled properties in mechanical metamaterials: A couple-stress theory perspective. Advanced Science 11, 2305113. doi:10.1002/advs.202305113.

Fang, H., Wu, H., Liu, Z., Zhang, Q., Xu, J., 2024. Evaluating dynamic models for rigid-foldable origami: unveiling intricate bistable dynamics of stacked-miura-origami structures as a case study. Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences 382, 20240014. doi:10.1098/rsta.2024.0014.

Feng, Y., Wang, M., Qiu, X., 2022. A simplified mechanical model of the crease in the flexible origami structures. International Journal of Solids and Structures 241, 111530. doi:10.1016/j.ijsolstr.2022.111530.

Filipov, E., Liu, K., Tachi, T., Schenk, M., Paulino, G., 2017. Bar and hinge models for scalable analysis of origami. International Journal of Solids and Structures 124, 26-45. doi:10.1016/J.IJSOLSTR.2017.05.028.

Filipov, E.T., Tachi, T., Paulino, G.H., 2015. Origami tubes assembled into stiff, yet reconfigurable structures and metamaterials. Proceedings of the National Academy of Sciences 112, 12321-12326. doi:10.1073/pnas.1509465112.

Gao, T., Bico, J., Roman, B., 2023. Pneumatic cells toward absolute Gaussian morphing. Science 381, 862-867. doi:10.1126/science.adi2997.

Gattas, J., You, Z., 2015. The behaviour of curved-crease foldcores under low-velocity impact loads. International Journal of Solids and Structures 53, 80-91. doi:10.1016/j.ijsolstr.2014.10.019.

Harris, J., McShane, G., 2021. Impact response of metallic stacked origami cellular materials. International Journal of Impact Engineering 147, 103730. doi:10.1016/j.ijimpeng.2020.103730.

Hassani, B., Hinton, E., 1998. A review of homogenization and topology optimization
i—homogenization theory for media with periodic structure. *Computers & Structures* 69, 707–717.
doi:10.1016/s0045-7949(98)00131-x.

Heimbs, S., 2013. Foldcore sandwich structures and their impact behaviour: An overview, in: *Dynamic Failure of Composite and Sandwich Structures*. Springer Netherlands, pp. 491–544.
doi:10.1007/978-94-007-5329-7_11.

Hollister, S.J., Kikuchi, N., 1992. A comparison of homogenization and standard mechanics analyses for periodic porous composites. *Computational mechanics* 10, 73–95. doi:10.1007/BF00369853.

Hu, Y., Liang, H., 2020. Folding simulation of rigid origami with lagrange multiplier method. *International Journal of Solids and Structures* 202, 552–561. doi:10.1016/j.ijsolstr.2020.06.016.

Hu, Y.C., Zhou, Y.X., Kwok, K.W., Sze, K.Y., 2021. Simulating flexible origami structures by finite element method. *International Journal of Mechanics and Materials in Design* 17, 801–829.
doi:10.1007/s10999-021-09538-w.

Jamalimehr, A., Mirzajanzadeh, M., Akbarzadeh, A., Pasini, D., 2022. Rigidly flat-foldable class of lockable origami-inspired metamaterials with topological stiff states. *Nature communications* 13, 1816. doi:10.1038/s41467-022-29484-1.

Jiang, B., Gao, W., Xie, Z., Pasini, D., Tan, H., 2025. Thin-walled tubular structures integrating origami patterns and tension-dominated bulkheads for enhanced energy absorption. *Thin-Walled Structures* 215, 113433. doi:10.1016/j.tws.2025.113433.

Jin, T., Li, L., Wang, T., Wang, G., Cai, J., Tian, Y., Zhang, Q., 2022. Origami-inspired soft actuators for stimulus perception and crawling robot applications. *IEEE Transactions on Robotics* 38, 748–764. doi:10.1109/TRO.2021.3096644.

Kalamkarov, A.L., Kolpakov, A.G., 1997. Analysis, design and optimization of composite structures. volume 1. Wiley New York. ISBN: 978-0-471-97189-4.

Kohn, R.V., Vogelius, M., 1984. A new model for thin plates with rapidly varying thickness. *International Journal of Solids and Structures* 20, 333–350. doi:10.1016/0020-7683(84)90044-1.

Lahiri, A., Pratapa, P.P., 2023. Folding-angle framework for structural modeling of rigid triangulated miura-ori lattices. *Journal of Mechanisms and Robotics* 15, 051004. doi:10.1115/1.4055742.

Lewinski, T., Telega, J., 1988. Asymptotic method of homogenization of two models of elastic shells. *Archiwum Mechaniki Stosowanej* 40, 705–723.

Li, J., Li, Q., Sun, T., Zhu, Z., Deng, Z., 2023. A general formulation for simulating the dynamic deployment of thick origami. *International Journal of Solids and Structures* 274, 112279. doi:10.1016/j.ijsolstr.2023.112279.

Li, M., Shen, L., Jing, L., Xu, S., Zheng, B., Lin, X., Yang, Y., Wang, Z., Chen, H., 2019a. Origami metawall: Mechanically controlled absorption and deflection of light. *Advanced science* 6, 1901434.
doi:10.1002/advs.201901434.

Li, S., Stampfli, J.J., Xu, H.J., Malkin, E., Diaz, E.V., Rus, D., Wood, R.J., 2019b. A vacuum-driven origami “magic-ball” soft gripper, in: 2019 International Conference on Robotics and Automation (ICRA), IEEE. pp. 7401–7408. doi:10.1109/ICRA.2019.8794068.

Liu, K., Paulino, G.H., 2017. Nonlinear mechanics of non-rigid origami: an efficient computational approach. *Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences* 473, 20170348. doi:10.1098/rspa.2017.0348.

Martinez, R.V., Fish, C.R., Chen, X., Whitesides, G.M., 2012. Elastomeric origami: programmable paper-elastomer composites as pneumatic actuators. *Advanced functional materials* 22, 1376–1384.
doi:10.1002/adfm.201102978.

Melancon, D., Gorissen, B., García-Mora, C.J., Hoberman, C., Bertoldi, K., 2021. Multi-stable inflatable origami structures at the metre scale. *Nature* 592, 545–550. doi:10.1038/
s41586-021-03407-4.

Mirzajanzadeh, M., Pasini, D., 2025. Reprogrammable curved-straight origami: Multimorphability and volumetric tunability. Science Advances 11, eadu4678. doi:10.1126/sciadv.adu4678.

Nassar, H., Lebée, A., Monasse, L., 2017. Curvature, metric and parametrization of origami tessellations: theory and application to the eggbox pattern. Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences 473, 20160705. doi:10.1098/rspa.2016.0705.

Nassar, H., Lebée, A., Werner, E., 2022. Strain compatibility and gradient elasticity in morphing origami metamaterials. Extreme Mechanics Letters 53, 101722. doi:10.1016/j.eml.2022.101722, arXiv:2207.08752.

Oudghiri-Idrissi, O., Guzina, B.B., 2022. Effective linear wave motion in periodic origami structures. Computer Methods in Applied Mechanics and Engineering 399, 115386. doi:10.1016/j.cma.2022.115386.

Overvelde, J.T.B., Weaver, J.C., Hoberman, C., Bertoldi, K., 2017. Rational design of reconfigurable prismatic architected materials. Nature 541, 347-352. doi:10.1038/nature20824.

Reddy, J., 1996. Mechanics of Laminated Composite Plates: Theory and Analysis. CRC-Press. ISBN: 0-8493-3101-3.

Reddy, J.N., 2006. Theory and analysis of elastic plates and shells. CRC press. doi:10.1201/9780849384165.

Schenk, M., Guest, S., 2011. Origami folding: A structural engineering approach, in: Origami 5. A K Peters/CRC Press, pp. 291-304. doi:10.1201/b10971-28.

Schenk, M., Guest, S.D., 2013. Geometry of miura-folded metamaterials. Proceedings of the National Academy of Sciences 110, 3276-3281. doi:10.1073/pnas.1217998110.

Sturm, R., Schatrow, P., Klett, Y., 2015. Multiscale modeling methods for analysis of failure modes in foldcore sandwich panels. Applied Composite Materials 22, 857-868. doi:10.1007/s10443-015-9440-9.

Turco, E., Barchiesi, E., dell'Isola, F., 2024. The long and winding road that leads to homogenisation of kresling origami. International Journal of Non-Linear Mechanics 163, 104756. doi:10.1016/j.ijnonlinmec.2024.104756.

Vasudevan, S.P., Pratapa, P.P., 2024. Homogenization of non-rigid origami metamaterials as kirchhoff-love plates. International Journal of Solids and Structures 300, 112929. doi:10.1016/j.ijsolstr.2024.112929.

Wang, C., Feng, L., Jasiuk, I., 2009. Scale and boundary conditions effects on the apparent elastic moduli of trabecular bone modeled as a periodic cellular solid. Journal of Biomechanical Engineering 131, 121008. doi:10.1115/1.4000192.

Wei, Z.Y., Guo, Z.V., Dudte, L., Liang, H.Y., Mahadevan, L., 2013. Geometric mechanics of periodic pleated origami. Physical review letters 110, 215501. doi:10.1103/PhysRevLett.110.215501.

Xue, W., Sun, Z., Ye, H., Liu, Q., Jian, B., Wang, Y., Fang, H., Ge, Q., 2024. Rigid-flexible coupled origami robots via multimaterial 3d printing. Smart Materials and Structures 33, 035004. doi:10.1088/1361-665X/ad212c.

Yang, Y., Read, H., Sbai, M., Zareei, A., Forte, A.E., Melancon, D., Bertoldi, K., 2024. Complex deformation in soft cylindrical structures via programmable sequential instabilities. Advanced Materials 36, 2406611. doi:10.1002/adma.202406611.

Yasuda, H., Miyazawa, Y., Charalampidis, E.G., Chong, C., Kevrekidis, P.G., Yang, J., 2019. Origami-based impact mitigation via rarefaction solitary wave creation. Science advances 5, eaau2835. doi:10.1126/sciadv.aau2835.

Zhang, J., Lu, G., Zhang, Y., You, Z., 2021. A study on ballistic performance of origami sandwich panels. International Journal of Impact Engineering 156, 103925. doi:10.1016/j.ijimpeng.2021.
103925.

Zhang, Q., Rudykh, S., 2024. Propagation of solitary waves in origami-inspired metamaterials. Journal of the Mechanics and Physics of Solids 187, 105626. doi:10.1016/j.jmps.2024.105626.

Zhang, X., Ma, J., Li, M., You, Z., Wang, X., Luo, Y., Ma, K., Chen, Y., 2022. Kirigami-based metastructures with programmable multistability. Proceedings of the National Academy of Sciences 119, e2117649119. doi:10.1073/pnas.2117649119.

Zheng, S., Yifeng, Z., Rong, L., Xiao, P., 2022. Vam-based reduced plate model for composite sandwich folded plate (CSFP) with v-shaped folded cores. Thin-Walled Structures 170, 108601. doi:10.1016/j.tws.2021.108601.