Programming shape using kirigami tessellations

TL;DR

This paper develops a novel inverse design framework for kirigami tessellations, enabling the creation of shape-changing sheets that conform to arbitrary 2D and 3D target shapes through optimized cut patterns.

Contribution

It introduces a constrained optimization approach to design kirigami patterns that deploy into specified shapes, extending capabilities to complex geometries and curved surfaces.

Findings

Validated designs with physical models confirming shape accuracy

Demonstrated deployment to various 2D target shapes

Extended framework to approximate 3D curved surfaces

Abstract

Kirigami tessellations, regular planar patterns formed by cutting flat, thin sheets, have attracted recent scientific interest for their rich geometries, surprising material properties and promise for technologies. Here we pose and solve the inverse problem of designing the number, size, and orientation of cuts that allows us to convert a closed, compact regular kirigami tessellation of the plane into a deployment that conforms approximately to any prescribed target shape in two and three dimensions. We do this by first identifying the constraints on the lengths and angles of generalized kirigami tessellations which guarantee that their reconfigured face geometries can be contracted from a non-trivial deployed shape to a novel planar cut pattern. We encode these conditions in a flexible constrained optimization framework which allows us to deform the geometry of periodic kirigami tesselations with three, four, and sixfold symmetry, among others, into generalized kirigami patterns that deploy to a wide variety of prescribed boundary target shapes. Physically fabricated models verify our inverse design approach and allow us to determine the tunable material response of the resulting structures. We then extend our framework to create generalized kirigami patterns that deploy to approximate curved surfaces in $\mathbb{R}^3$. Altogether, this work illustrates a novel framework for designing complex, shape-changing sheets from simple cuts showing the power of kirigami tessellations as flexible mechanical metamaterials.

Equations

$$\Vert {\mathbf{p}}_i-{\tilde{\mathbf{p}}}_i{\Vert }^2=0,$$

$$a^2-b^2=0.$$

$${\theta }_1+{\theta }_2+{\theta }_3+{\theta }_4=2\pi ,$$

$$⟨(\mathbf{b}-\mathbf{a})\times (\mathbf{c}-\mathbf{a}),\vec{n}⟩\ge 0,$$

$$\frac{1}{M}\sum _{i=1}^M\left(\sum _j({\alpha }_{i_j}-{\beta }_{i_j}{)}^2+\sum _k(a_{i_k}-b_{i_k}{)}^2\right)$$

$$E({\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_N)=\frac{1}{N_s}\sum _{i,j}{\left(\frac{\Vert {\mathbf{x}}_i-{\mathbf{x}}_j\Vert -l_{ij}}{l_{ij}}\right)}^2+\lambda \frac{1}{N_c}\sum _i{\theta }_i^2,$$

$$\Vert {\mathbf{x}}_i-{\tilde{\mathbf{x}}}_i{\Vert }^2=0,$$

$$⟨(\mathbf{b}-\mathbf{a})\times (\mathbf{c}-\mathbf{a}),(\mathbf{c}-\mathbf{a})\times (\mathbf{d}-\mathbf{a})⟩\ge 0,$$

$$\text{Volume}(F)=0.$$

$$⟨(\mathbf{b}-\mathbf{a})\times (\mathbf{c}-\mathbf{a}),\mathbf{d}-\mathbf{a}⟩=0,$$

Full text

Programming shape using kirigami tessellations

Gary P. T. Choi1, Levi H. Dudte1, L. Mahadevan1,2,3∗

1John A. Paulson School of Engineering and Applied Sciences, Harvard University

2Departments of Physics, and Organismic and Evolutionary Biology, Harvard University

3 Kavli Institute for Nanobio Science and Technology, Harvard University, Cambridge, MA 02138, USA

∗To whom correspondence should be addressed; E-mail: lmahadev@g.harvard.edu

Kirigami tessellations, regular planar patterns formed by cutting flat, thin sheets, have attracted recent scientific interest for their rich geometries, surprising material properties and promise for technologies. Here we pose and solve the inverse problem of designing the number, size, and orientation of cuts that allows us to convert a closed, compact regular kirigami tessellation of the plane into a deployment that conforms approximately to any prescribed target shape in two and three dimensions. We do this by first identifying the constraints on the lengths and angles of generalized kirigami tessellations which guarantee that their reconfigured face geometries can be contracted from a non-trivial deployed shape to a novel planar cut pattern. We encode these conditions in a flexible constrained optimization framework which allows us to deform the geometry of periodic kirigami tesselations with three, four, and sixfold symmetry, among others, into generalized kirigami patterns that deploy to a wide variety of prescribed boundary target shapes. Physically fabricated models verify our inverse design approach and allow us to determine the tunable material response of the resulting structures. We then extend our framework to create generalized kirigami patterns that deploy to approximate curved surfaces in $\mathbb{R}^{3}$. Altogether, this work illustrates a novel framework for designing complex, shape-changing sheets from simple cuts showing the power of kirigami tessellations as flexible mechanical metamaterials.

Kirigami is the creative art of paper cutting and folding that originated in Japan. Kirigami tessellations, which are regular planar patterns formed by cuts, have recently emerged as paradigms of mechanical metamaterials. Various studies have focused on their geometry and kinematics  (?, ?, ?), the mechanics of their deployment  (?, ?, ?, ?, ?, ?), the mathematics of their constructions  (?), and as materials with a range of unusual properties such as topological insulators  (?, ?), auxetic structures  (?, ?), super-stretchable sheets  (?, ?, ?, ?, ?) etc. However, almost without exception, these studies study the forward problem of understanding the behavior given the topology and geometry of the tessellations. From both a mathematical and a technological perspective, a very natural question that arises is this: how can one design the geometry of the cuts in a given planar structure, i.e. tessellate a closed, compact subset of the plane, so that it can be deployed into a prescribed final shape in two or three dimensions? Here, we answer this question using an inverse approach couched as a constrained optimization problem, that has some algorithmic similarities to our recent solution of the inverse problem for origami tessellations for programming curvature  (?).

We start from the well known observation that the only regular polygons that can tile the plane are triangle, square and hexagon. While any of these simple tilings and even more complex ones  (?, ?) can be the basis for a deployable kirigami pattern, to simplify our discussion, we present our framework using the quad kirigami pattern (see Section S1 in Supplementary Information - SI for the others). In Fig. 1a, we show the quadrilateral tessellation of the plane in its compact and deployed state, with the cuts along the edges of the quads designed to allow for rotational in-plane deployment about hinges. The deployment associated with thus 1-degree of freedom mechanism yields a continuous family of self-similar shapes that is geometrically self-limiting eventually. The basic unit cell underlying this pattern is also shown in Fig. 1a in both the undeployed and deployed states and clearly shows the sides that must be identified together and the hinges about which the quads rotate, both of which provide constraints on the class of admissible solutions.

Our inverse design problem then becomes this: how should the unit cell be modulated in space to approximate a given planar shape in its deployed state and still be able to tile the plane compactly when undeployed. In Fig. 1b, we outline this procedure. To find a pattern that can deploy to approximate a prescribed planar shape, it is natural to perform the search for the admissible results in the deployed space. Therefore, we start with some approximation to the given deployed configuration that matches the prescribed boundary curve that might be obtained using, for example, a conformal/quasi-conformal map  (?, ?) (although we could use any approximation that preserves the number and connectedness of the quads). In general the deployed configuration will not correspond to a deployment of a kirigami pattern as the associated length and angle constraints at every interior node are not satisfied. Imposing these interior and boundary conditions leads to a constrained optimization problem whose solution naturally leads to an admissible deployed kirigami pattern whose compact tiling follows via a simple, albeit non-affine contraction.

To ensure that for a given boundary curve $\partial S$ we have a valid deployed configuration, we formulate the constraint as the boundary shape matching constraints. The constraints force all nodes on the boundary of the deployed configuration to lie exactly on $\partial S$. Mathematically, for every boundary node $\mathbf{p}_{i}$, we should have

[TABLE]

where $\widetilde{\mathbf{p}}_{i}$ is the projection of $\mathbf{p}_{i}$ onto $\partial S$ and $\|\cdot\|$ is the Euclidean 2-norm. In addition to matching the target boundary shape in the deployed configuration, we can also control the boundary shape of the kirigami tessellations in the undeployed state, although here we neglect this additional condition (see Section S2 in SI for details of how to incorporate this as well).

For a valid deployed configuration, we must also be able to contract (undeploy) the configuration into a generalized kirigami pattern that closes consistently along the cuts without any mismatch or overlap in lengths and angles. As illustrated in Fig. 1b, for a valid deployed configuration of a generalized kirigami pattern, the following contractibility constraints should be satisfied: (i) for every pair of edges with edge lengths $a,b$ in the deployed space that correspond to the same cut,

[TABLE]

(ii) For every set of four angles in the deployed space that correspond to an interior node, their sum should be $2\pi$:

[TABLE]

where $\theta_{i}$ are angles in the deployed space as illustrated in Fig. 1a. The analogous formulation for other kirigami tessellations are discussed in Section S2 in SI.

While the constraints described above ensure a consistency between corresponding edges and angles, they do not prevent the faces from overlapping. To enforce this, we use the following non-overlap constraints at every angle between two adjacent faces:

[TABLE]

where $\mathbf{a}$ and $\mathbf{b}$ are two nodes of a face, $\mathbf{a}$ and $\mathbf{c}$ are two nodes of another face, $(\mathbf{b},\mathbf{a},\mathbf{c})$ form a positive (right-hand ordered) angle between the two faces, and $\vec{n}=(0,0,1)$ is the outward unit normal.

Solving a constrained optimization problem involving all the above constraints is sufficient to guarantee a valid deployed configuration of a generalized kirigami pattern that approximates the prescribed shape. However, the solution is likely to be very rough with large gradients in the shapes of the quads. To get a smooth kirigami tessellation, we therefore minimize the following objective function:

[TABLE]

where $\alpha_{i_{j}},\beta_{i_{j}}$ are a pair of corresponding angles in two adjacent cells and $a_{i_{k}},b_{i_{k}}$ are corresponding edge lengths in two adjacent cells, and $M$ is the total number of pairs of adjacent cells, subject to the constraints (1),(2),(3),(4). We solve the problem numerically using MATLAB’s built-in optimization routine fmincon (see Section S3 in SI for details). We note that our optimization problem is underconstrained so that there are multiple admissible deployable kirigami patterns (see Section S2 and S4 in SI for details). Once we find an admissible deployed kirigami pattern, we can contract this into its compact form by rotating the faces contracting the entire structure.

To illustrate the effectiveness of our inverse design approach, we show how to circle the square and how to make an egg shape from a square via generalized kirigami patterns. As shown in Fig. 2a, by introducing novel generalized kirigami patterns on a square, the deployed configurations can effectively match either a circle or an egg.

Our method can also generate novel generalized kirigami patterns that, when deployed, approximate boundary shapes with different curvatures. Fig. 2b shows two generalized kirigami patterns; one deploys to approximate a boundary shape with mixed curvature, and the other deploys to approximate a rectangle (see Section S4 in SI for more generalized kirigami patterns with other base tessellations). These generalized kirigami patterns with different base tessellations exhibit different behaviors in terms of their porosity and magnifying property (see Section S5 in SI for details).

There is naturally an accuracy-effort trade-off in approximating a prescribed shape using generalized kirigami tessellations: a larger number of smaller tiles allows for improved approximation at the cost of more tiles. Fig. 2c shows several generalized kirigami patterns of circling the square with different resolution; with more tiles, the boundary of the deployed pattern gets closer to a perfect circle (see Section S4 in SI for multiresolution results for other patterns). To quantitatively assess the accuracy of the patterns for the approximation, we define the boundary layer area of a generalized kirigami pattern by the total area of the gaps between the target boundary shape and the boundary of the deployed configuration. From the log-log plot in Fig. 2c, it can be observed that as the boundary layer area (denoted by $A$) decreases as the number of tiles (denoted by $n^{2}$) increases. More specifically, the slope of the least-square line is approximately $-1/2$, which indicates that $A\propto(n^{2})^{-1/2}=n^{-1}$. To further explain this, we approximate every boundary gap by a triangle and measure the change in the average triangle base length $\tilde{l}$ and average triangle height $\tilde{h}$ for different resolutions. We observe that $\tilde{l}\propto n^{-1}$ and $\tilde{h}\propto n^{-1}$, and hence the average area of the triangles $\tilde{a}\propto n^{-2}$. As the number of boundary gap triangles is approximately $4n$, we have $A\approx 4n\tilde{a}\propto n^{-1}$.

Our inverse design approach guarantees that the end-points, i.e. the deployed and undeployed states satisfy the contractibility and boundary constraints, but are agnostic to whether there is a one-parameter family of zero-energy solutions linking them. To see if this is indeed true, we simulate the deployment process from the initial state to the deployed state by considering a simple physical model of the kirigami tessellation using linear springs along the edges and diagonals of the quads, and simple torsional springs at the nodal hinges to model the ligaments that hold the structure together in real systems. Then, the mechanical energy of the system is given by

[TABLE]

where $\mathbf{x}_{i}$ are the locations of the nodes, $\theta_{i}$ are the angles between every pair of edges created under the cuts, $l_{ij}$ are the rest lengths of the extensional springs, $N_{s}$ is the total number of extensional springs, $N_{c}$ is the total number of torsional springs, and $\lambda$ is the ratio of the torsional spring constant to the extensional spring constant. A larger $\lambda$ corresponds to a thicker ligament, which has a larger tendency to close up. By iteratively moving the boundary nodes to the target boundary shape and solving for the intermediate deployed configurations, we obtain a continuous deployment simulation. Fig. 3a shows the energetics of the deployment simulations with different $\lambda$: if $\lambda\to 0$, we see the presence of bistability, while if $\lambda\not\to 0$, monostability or multistability can be observed. To verify these results, we fabricated a physical model by laser cutting sheets of super-stretchable abrasion-resistant natural rubber. Fig. 3b shows the deployment snapshots of a fabricated model with monostability, showing that the computational design and the real deployment have similar behaviors (see Section S4 in SI for another example).

While our inverse design approach was originally focused on approximating planar shapes, it can be extended to fit surfaces in three dimensions, i.e. the pattern space is in $\mathbb{R}^{2}$ while the deployed space is in $\mathbb{R}^{3}$. To fit a surface $S$ in $\mathbb{R}^{3}$ instead of a given planar boundary shape, we replace the boundary shape matching constraints (1) by the surface matching constraints below. For every node $\mathbf{x}_{i}$ in the deployed configuration, we enforce

[TABLE]

where $\widetilde{\mathbf{x}}_{i}$ is the projection of $\mathbf{x}_{i}$ onto $S$ and $\|\cdot\|$ is the Euclidean 2-norm. The extra constraints for controlling the boundary shape of the undeployed configuration in the planar case can be directly extended to the surface fitting problem.

One can easily note that the contractibility constraints for surface fitting are the same as (2) and (3). For the non-overlap constraints, to prevent adjacent faces in the deployed configuration from overlapping or intersecting between each other, we enforce the following inequality constraints for every two adjacent faces in the deployed configuration:

[TABLE]

where $\mathbf{a},\mathbf{b}$ are two nodes of a face, $\mathbf{c},\mathbf{a},\mathbf{d}$ are three nodes of another face, $(\mathbf{b},\mathbf{a},\mathbf{c})$ form a positive (right-hand ordered) angle between the two faces, and $(\mathbf{c},\mathbf{a},\mathbf{d})$ also form a positive (right-hand ordered) angle. The idea is to replace the unit normal $\vec{n}$ in (4) by the normal computed using one of the two faces.

Note that in the case of deployed configurations for planar shape matching, all faces are automatically planar. This property does not necessarily hold for surface fitting. In addition to the above constraints, we need to enforce the following planarity constraints for each face $F$:

[TABLE]

More explicitly, for quad tessellations the constraint becomes

[TABLE]

where $\mathbf{a},\mathbf{b},\mathbf{c},\mathbf{d}$ are the four points of the quad $F$.

Finally, both the objective function (5) for the planar case and the contraction process can be easily extended for surface fitting. Hence, we can obtain a valid generalized kirigami pattern that deploys to approximate a prescribed surface by solving a constrained optimization problem analogous to the planar one using fmincon in MATLAB, with (7),(2),(3),(8),(9) to be satisfied.

Fig. 4 shows several generalized kirigami patterns that deploy to fit surfaces with different curvature properties. Again, we can impose extra constraints to further achieve different effects. For instance, we can obtain a rectangular quad pattern that deploys to fit a hyperbolic paraboloid.

It is noteworthy that surfaces in $\mathbb{R}^{3}$ have non-zero curvature in general, while all faces of the deployed configurations of the generalized kirigami patterns must have zero curvature. This suggests that the curvature of the surfaces resides in the holes of the deployed configurations of generalized kirigami patterns. To verify this, we consider fitting every hole in the deployed configurations of generalized quad kirigami patterns by a bicubic Bézier surface, and compute the mean curvature and the Gaussian curvature of it (see Section S6 in SI for the details). From the curvature plots in Fig. 4, it can be observed that the curvatures of the underlying surfaces indeed reside at the holes of the deployed generalized kirigami patterns, especially at the highly curved regions of the underlying surfaces.

Our inverse design approach allows us to create generalized kirigami patterns that approximate a large variety of planar shapes when deployed. Deployment simulations and fabricated models have shown that the generalized kirigami patterns achieve monostability with thick hinges, and by reducing the size of the hinges we can achieve bistability. Our method can also be extended for fitting surfaces in 3D. By establishing a framework for designing generalized kirigami patterns, our work has set the stage for a better understanding of the interplay between geometry, topology and mechanics.

Acknowledgment This work was supported in part by the Croucher Foundation (to G.P.T.C.), National Science Foundation Grant DMR 14-20570 (to L.H.D. and L.M.) and DMREF 15-33985 (to L.H.D. and L.M.).

Competing Interests We have filed a patent on our algorithms for kirigami design.

Data Availability Additional information is available in the Supplementary Information. Correspondence and requests for materials should be addressed to L.M..

References and Notes

Grima, J. N., & Evans, K. E.. Auxetic behavior from rotating squares. J. Mater. Sci. Lett. 19, 1563-1565 (2000).

• 2.

Grima, J. N., Alderson, A., & Evans, K. E.. Negative Poisson’s ratios from rotating rectangles. Comp. Methods Sci. Technol. 10, 137-145 (2004).

• 3.

Grima, J. N., Alderson, A., & Evans, K. E.. Auxetic behaviour from rotating rigid units. Phys. Status Solidi B 242, 561-575 (2005).

• 4.

Rafsanjani, A., & Pasini, D.. Bistable auxetic mechanical metamaterials inspired by ancient geometric motifs. Extreme Mech. Lett. 9, 291-296 (2016).

• 5.

Sussman, D. M., Cho, Y., Castle, T., Gong, X., Jung, E., Yang, S., & Kamien, R. D.. Algorithmic lattice kirigami: a route to pluripotent materials. Proc. Natl Acad. Sci. USA 112, 7449-7453 (2013).

• 6.

Blees, M. K., Barnard, A. W., Rose, P. A., Roberts, S. P., McGill, K. L., Huang, P. Y., Ruyack, A. R., Kevek, J. W., Kobrin, B., Muller, D. A. & McEuen, P.L.. Graphene kirigami. Nature 524, 204-207 (2015).

• 7.

Zhang, Y., Yan, Z., Nan, K., Xiao, D., Liu, Y., Luan, H., Fu, H., Wang, X., Yang, Q., Wang, J. & Ren, W.. A mechanically driven form of Kirigami as a route to 3D mesostructures in micro/nanomembranes. Proc. Natl Acad. Sci. USA 112, 11757-11764 (2015).

• 8.

Shyu, T.C., Damasceno, P.F., Dodd, P.M., Lamoureux, A., Xu, L., Shlian, M., Shtein, M., Glotzer, S.C. & Kotov, N.A.. A kirigami approach to engineering elasticity in nanocomposites through patterned defects. Nat. Mater. 14, 785 (2015).

• 9.

Rafsanjani, A. & Bertoldi, K.. Buckling-induced kirigami. Phys. Rev. Lett. 118, 084301 (2017).

• 10.

Konaković, M., Crane, K., Deng, B., Bouaziz, S., Piker, D., & Pauly, M.. Beyond developable: computational design and fabrication with auxetic materials. ACM Trans. Graph. 35, 89 (2016).

• 11.

Kane, C. L., & Lubensky, T. C.. Topological boundary modes in isostatic lattices. Nat. Phys. 10, 39 (2014).

• 12.

Sun, K., Souslov, A., Mao, X., & Lubensky, T. C.. Surface phonons, elastic response, and conformal invariance in twisted kagome lattices. Proc. Natl Acad. Sci. USA 109, 12369-12374 (2012).

• 13.

Tang, Y., & Yin, J.. Design of cut unit geometry in hierarchical kirigami-based auxetic metamaterials for high stretchability and compressibility. Extreme Mech. Lett. 12, 77-85 (2017).

• 14.

Gatt, R., Mizzi, L., Azzopardi, J. I., Azzopardi, K. M., Attard, D., Casha, A., Briffa, J. & Grima, J. N.. Hierarchical auxetic mechanical metamaterials. Sci. Rep. 5, 8395 (2015).

• 15.

Kolken, H. M., & Zadpoor, A. A.. Auxetic mechanical metamaterials. RSC Adv. 7, 5111-5129 (2017).

• 16.

Mitschke, H., Robins, V., Mecke, K., & Schröder-Turk, G. E.. Finite auxetic deformations of plane tessellations. Proc. R. Soc. Lond. A 469, 20120465 (2013).

• 17.

Shan, S., Kang, S. H., Zhao, Z., Fang, L., & Bertoldi, K.. Design of planar isotropic negative Poisson’s ratio structures. Extreme Mech. Lett. 4, 96-102 (2015).

• 18.

Isobe, M., & Okumura, K.. Initial rigid response and softening transition of highly stretchable kirigami sheet materials. Sci. Rep. 6, 24758 (2016).

• 19.

Neville, R. M., Scarpa, F. & Pirrera, A.. Shape morphing Kirigami mechanical metamaterials. Sci. Rep. 6, 31067 (2016).

• 20.

Dudte, L. H., Vouga, E., Tachi, T., & Mahadevan, L.. Programming curvature using origami tessellations. Nat. Mater. 15, 583 (2016).

• 21.

Grünbaum, B., & Shephard, G. C.. Tilings and patterns. Freeman (1987).

• 22.

Choi, G. P.-T., & Lui, L. M.. A linear formulation for disk conformal parameterization of simply-connected open surfaces. Adv. Comput. Math. 44, 87-114 (2018).

• 23.

Meng, T. W., Choi, G. P.-T., & Lui, L. M.. TEMPO: feature-endowed Teichmüller extremal mappings of point clouds. SIAM J. Imaging Sci. 9, 1922-1962 (2016).