Research ArticleENGINEERING

Origami-based impact mitigation via rarefaction solitary wave creation

See allHide authors and affiliations

Science Advances  24 May 2019:
Vol. 5, no. 5, eaau2835
DOI: 10.1126/sciadv.aau2835


The principles underlying the art of origami paper folding can be applied to design sophisticated metamaterials with unique mechanical properties. By exploiting the flat crease patterns that determine the dynamic folding and unfolding motion of origami, we are able to design an origami-based metamaterial that can form rarefaction solitary waves. Our analytical, numerical, and experimental results demonstrate that this rarefaction solitary wave overtakes initial compressive strain waves, thereby causing the latter part of the origami structure to feel tension first instead of compression under impact. This counterintuitive dynamic mechanism can be used to create a highly efficient—yet reusable—impact mitigating system without relying on material damping, plasticity, or fracture.


Mechanical metamaterials offer a new dimension in achieving nonconventional and tailored mechanical properties through architecture (13). In particular, the manipulation of wave propagation in mechanical metamaterials is a topic of intense research for various engineering applications, e.g., waveguiding, vibration filtering, subwavelength imaging, and impact mitigation (47). The ability to achieve this desirable mechanical performance often relies on the platform in which we construct mechanical metamaterials. The choice of platform can vary from periodically arranged microlattice/macrolattice structures to self-assembling particles to three-dimensional (3D) printed soft/hard architected materials (811).

Recent studies have shown that origami can serve as an ideal playground to realize highly versatile and tunable mechanical metamaterials (1216). For example, by introducing crease lines into flat surface materials, one can construct origami-based structures, which can offer enhanced stiffness (17), negative Poisson’s ratio (18), and multistability (1921). Given the scale-free nature of origami, this kind of design framework can be used in a wide range of scales. For example, the concept of origami has been adapted to a diverse set of design principles, including robotics (22), reconfigurable structures (23), and self-folding actuated by living cells (24). While these studies focused mainly on origami’s static or quasi-static behavior, the analysis of the dynamics of origami-based structures is a natural next step to investigate. However, the connection between the origami crease pattern and the dynamic folding/unfolding behavior of origami itself has been relatively unexplored (25, 26). In particular, very few experimental studies have been reported (27, 28).

In the present study, we explore unique wave dynamics in a mechanical metamaterial that is composed of volumetric origami structures. In particular, each volumetric origami structure is a triangulated cylindrical origami (TCO) (Fig. 1A for the folding motion and Fig. 1B for the flat crease pattern before it is assembled into a unit cell) (21). This TCO unit cell is analogous to a post-buckled shape of a cylinder under simultaneous axial and twisting loading (29, 30). Thus, the folding motion of the TCO is characterized by coupling between axial and rotational motions, as shown in Fig. 1A. Recent studies on the TCO structure have shown their versatile mechanical properties, such as its tailorable stability (20), zero-stiffness mode (31), and strain-softening/hardening behavior (21).

Fig. 1 Geometry of the TCO prototypes.

(A) Folding motion of the TCO is shown in sequence. (B) The flat sheet with crease patterns (upper left) is composed of mountain crease lines (red), valley crease lines (blue), and the adhesive area (shaded area). The photograph shows corresponding laser-cut paper sheets (lower right). (C) Actual prototype of the origami-based metamaterial and its unit cell (lower right inset). (D) The origami-based metamaterial generates the rarefaction solitary wave despite the application of compressive impact. The system is composed of the TCO unit cells (lower right). To connect the neighboring unit cells, we use the interfacial polygonal cross-section with markers at vertices (lower left). Photo credit: H.Y. and Y.M., University of Washington.

Our origami structure consists of 20 TCO unit cells, which are fabricated by using paper sheets cut by a laser cutting machine (Fig. 1, B and C; see also Materials and Methods for details). The crease patterns are carefully designed to make the TCO cells exhibit effective strain-softening behavior, as in (21). This softening behavior in origami is in sharp contrast to conventional nonlinear metamaterials that mainly exploit strain-hardening behavior, such as granular crystals (8, 11, 32), that have been serving as a popular test bed for demonstrating classical dynamics of nonlinear mass-spring systems [e.g., the celebrated Fermi-Pasta-Ulam-Tsingou chain (33)].

Recently, it has been theoretically and numerically predicted that a 1D discrete system with strain-softening interactions can support the propagation of a solitary wave in the form of a rarefaction pulse (8, 34, 35). In this 1D system, the application of a compressive impact generates a tensile solitary wave, propagating ahead of the initial compressive strain wave and thereby making the latter part of the medium feel tension first instead of compression (see Fig. 1D for the conceptual illustration). While the feasibility of verifying this counterintuitive dynamics has been discussed in different settings (35, 36), the experimental demonstration of these rarefaction solitary waves remains elusive. This is mainly due to the challenges in fabricating an effective strain-softening medium that exhibits minimal dissipation for waveguiding purposes.

Here, we report the first experimental observation of a mechanical rarefaction solitary wave in the strain-softening setting of TCO-based metamaterials. This rarefaction wave overtakes the leading compressive wave generated by a compressive impact. The experimentally observed rarefaction pulse agrees quantitatively with numerical simulations of the derived model and agrees qualitatively with an asymptotic analysis based on the celebrated Korteweg–de Vries (KdV) equation (37, 38).


We start by examining the geometry of the TCO unit cell. The TCO structure has axial and rotational motions that are coupled to each other. We characterize this folding motion by using the axial displacement (u) and rotational angle (φ), which are defined with respect to the initial height (h0) and angle (θ0) (Fig. 2 , A and B). To analyze the kinematics of the origami, we construct a two–degree of freedom (2DOF) mathematical model, which approximates the folding behavior of the TCO into linear spring motions along the crease lines (21). In this model, we use two different spring constants (Ka and Kb) for the shorter (AaB¯ in Fig. 2A) and longer crease lines (AbB¯), respectively. These two linear springs are repeated along the Np-sided polygon (denoted by the shaded area with the radius R of the circle circumscribing the polygon). In addition, the side crease (AaAb¯) is modeled as a torsional spring to capture the folding of the side facets with the angle of ψ (Fig. 2A). The torsional spring constant for this crease line is denoted by Kψ. These three spring constants are determined empirically by conducting compression tests on the fabricated TCO cells (see Materials and Methods for details and movie S1 for the fabrication and folding motion of the TCO).

Fig. 2 Folding motions of the TCO with strain-softening behavior.

(A) The axial displacement (u) is defined with respect to the initial height (h0) of the TCO. (B) Top-down view shows the rotational angle (ϕ) defined with respect to the initial angle (θ0). (C) Axial force (F normalized by the spring constant Ka and h0) versus displacement. The dashed red curve with the colored area represents the experimental value with the SD. The solid blue curve denotes the 2DOF linear spring model. The inset shows the variations of the stiffness as a function of the TCO displacement.

We analyze the elastic potential energy change of the TCO unit cell as a function of u and φ as followsU(u,φ)=12NpKa(aa(0))2+12NpKb(bb(0))2+12(2Np)Kψ(ψψ(0))2(1)where a = a(u, φ) and b = b(u, φ) are the length of the crease lines (AaB¯ and AbB¯), respectively, and the superscript (0) denotes the initial states. By using this energy expression and applying the minimum potential energy principle (21), we obtain the 2DOF model that indicates how the deformation of the TCO cell takes place by coupling axial and rotational motions (see the “Equations of motion” section in the Supplementary Materials for details). Figure 2C shows the analytical prediction (solid curve) obtained from this 2DOF model, along with the experimental force and displacement relationship (dashed curve). See Materials and Methods for the detailed testing and table S1 for the parameters used in experiments and analytics. We observe the stiffness of the structure decreases as it is compressed, while the trend is opposite under tension (inset of Fig. 2C). We will later discuss how this strain-softening behavior contributes to the unique wave dynamics in this TCO-based system.


With a firm understanding of the unit cell at hand, we are ready to investigate the wave dynamics in a chain of TCO cells (Fig. 3A). For experiments, we prototype 20 identical TCO unit cells by choosing h0 = 35 mm, θ0 = 70, R = 36 mm, and Np = 6. To ensure uniform and repeatable folding behavior in these units, we apply the preconditioning process to each cell and verify the uniformity by measuring its compressive motions (see Materials and Methods for details). The left end of the chain (unit number n = 1) is connected to a shaker through the customized attachment with a sleeve bearing, which transfers the shaker impact to the cell while allowing its free rotational motion (see the upper inset of Fig. 3A). The TCO cell positioned in the right end of the chain (n = 20) is fixed to the rigid wall. To measure the dynamic folding/unfolding motion of each unit cell, we use the digital image correlation (DIC) technique by using three pairs of action cameras (GoPro HERO4 Black) whose maximum frame rate is 240 frames per second (see the lower inset of Fig. 3A).

Fig. 3 Experimental setup and DIC analysis results.

(A) The shaker is attached to the leftmost unit cell through the sleeve bearing (upper left inset). The folding motion of each unit cell is captured by six action cameras (lower inset). For DIC analysis, the fluorescent green markers are used. (B) Snapshots of the experiment at t = 0, 0.06, 0.11, and 0.14 s. Images from the camera are shown in the left column, where the red (blue) arrows represent the compressive (tensile) velocity vector of the polygon in the axial direction. 3D reconstruction of the TCO chain (right column). The deformation is scaled 2.5 times larger than the original deformation for visual clarity. The gray arrows indicate the propagation of the rarefaction solitary wave. Photo credit: H.Y. and Y.M., University of Washington.

The digital images in Fig. 3B show the snapshots of the first eight TCO unit cells at four different time frames captured by the first pair of the action cameras (see movie S2). In this figure, the length of the colored arrows represents the speed of the polygon in the axial direction, and red (blue) color denotes rightward (leftward) velocity. For the sake of visualization, the 3D images are reconstructed based on the experimental data as shown in the right column of Fig. 3B, where red (blue) color indicates compressive (tensile) strain (see movie S3). Here, the strain value of the nth cell is defined by (unun + 1)/h0, where un (un + 1) is the axial displacement of its left (right) polygon such that compressive strains take positive signs for convenience. From the experimental results, we observe that, initially, the first unit shows the large-amplitude compression due to the excitation by the shaker. However, this compressive motion decays quickly without being robustly transmitted along the chain, whereas the noticeable tensile motion is evolved instead (see the gray arrows in Fig. 3B and also movie S2). Thus, it appears that a tensile wave has propagated despite the application of a compressive force.

To conduct a more thorough analysis of this counterintuitive wave dynamics, we plot the measured strains in time and space domains (Fig. 4A). We observe evidently that the last TCO cell (n = 20) experiences a tensile strain (see the black arrow) despite the application of compressive impact to the system. This tensile strain is due to the formation and propagation of a rarefaction solitary wave (blue valley as indicated by the green arrow). Rarefaction (i.e., tensile) solitary waves are the conical type of solitary wave in strain-softening systems (34), such as the TCO chain considered herein. One would not expect to find compressive solitary waves in our system because they are observed typically in strain-hardening systems, such as granular crystals (8).

Fig. 4 Wave form analysis.

(A) Space-time evolution of the experimentally measured strain wave propagation in the origami-based system. The black arrow indicates the rarefaction solitary wave, and the green one shows the direction of the propagation. (B) Numerical simulation results show a qualitative agreement with the experimental data. The black arrow indicates the leading compressive wave in front of the rarefaction wave. (C) Amplitude change of the rarefaction solitary wave. The experimental data are fitted by the KdV solution (black curve) to obtain the damping coefficient for numerical and analytical analysis. The error bar represents the SD calculated from five measurements. Simulation results are shown in blue dots, which are fitted to the blue dashed curve. (D) The amplitude of the leading compression is analyzed. The dashed curves are obtained from the exponential fit to the experimental and numerical data. The inset shows the exponential decay of the compressive strain. The shapes of the rarefaction solitary wave (E) at t = 0.10 s and (F) t = 0.15 s are shown.

To better understand the above experimental results, we formulate the equations of motion for the entire chain based on the 2DOF model of the TCO cells and solve the resulting equations numerically by using the fourth-order Runge-Kutta method (see the “Equations of motion” section in the Supplementary Materials for details). The simulation results are shown in Fig. 4B, which demonstrates excellent agreement with the experimental results. The dissimilar trend of wave propagation between the tensile and compressive waves in our system is worth noting (compare the blue and red peaks indicated by black arrows in Fig. 4, A and B). To examine the difference between these two types of waves, we extract and compare the strains of the leading tensile and compressive waves in Fig. 4, C and D, respectively. We find that the compressive component decays more drastically than the tensile counterpart. The compressive waves show an order-of-magnitude reduction in amplitude, as verified by the exponential decay of the compressive strain in the inset of Fig. 4D.

The strong contrast of wave attenuation trends stems from the intrinsic nature of the two different waves formed in the tensile and compressive realms. In the tensile domain, our TCO system generates rarefaction solitary waves in a robust form, while in the compressive region, it forms oscillatory waves that disperse energy into space and time (see oscillations in Fig. 4, A and B). As will be verified later, the rarefaction solitary waves are supersonic and thus nonlinear, while the oscillatory ones are linear and limited by the sound speed. These dissimilar mechanisms are attributed to the strain-softening nature of the TCO system in compression, where the cell stiffness decreases (increases) as a further compression (tension) is imposed (see the inset of Fig. 2C). In nonlinear mass-spring systems, this asymmetric behavior implies that the system will favor wave localization in tension and dispersion in compression (34).

Although general strain-softening systems can also be used for impact mitigation purposes, they typically require hundreds of cells to generate rarefaction solitary waves based on previous numerical studies (3436). In the present experimental work, however, we find that our TCO architecture can form rarefaction waves even within 10 TCO cells by leveraging the tailorable strain-softening nature of the TCO and its unique wave-coupling motions. This manifests the efficacy of the TCO-based metamaterial in mitigating the compressive impact.

The robust propagation of the leading tensile waves can also be modeled by theoretical analysis. To this end, we first postulate unn based on the eigenmode of the single TCO unit cell. Note that this linear relationship holds approximately in the full range of dynamic strains considered here (see the “Derivation of single-component model” section in the Supplementary Materials for an a posteriori validation of this approximation). On the basis of this, the equations of motion can be reduced to a single-component model. Note that the single-component model is still nonlinear by virtue of the force-displacement and torque-angle relationships. Then, we take the continuum limit of this single-component equation in the infinite TCO chain to derive the well-known KdV equation (see the “Continuum limit” section in the Supplementary Materials). This equation has a closed-form rarefaction solitary wave solution, thereby yielding an analytical approximation of the wave in the TCO lattice. Note that this analytical approach captures the effect of damping (see the gradual decrease of the KdV curve in Fig. 4C) based on the empirical dashpot model (see the “Equation of motion with damping effect” and “Continuum limit” sections in the Supplementary Materials).

By using the experimental, numerical, and analytical data, we now investigate the evolution of the strain waves at t = 0.10 and 0.15 s, as shown in Fig. 4 (E and F). In these plots, we already witness a substantially attenuated form of the leading compressive wave. On the contrary, the tensile one is propagating more dominantly and robustly. Although the analytical waveform focuses only on the tensile components, we observe a qualitative agreement among analytical, numerical, and experimental results. One interesting finding here is that the maximum compressive strain is ahead of the rarefaction solitary wave at t = 0.10 s, but the peak is shifted to the rear part of the solitary wave at t = 0.15 s. This indicates that the rarefaction solitary wave propagates faster than the other compressive wave packets and eventually overtakes the initial compressive strain wave. While this overtaking behavior has been previously studied numerically (34, 35), this is the first experimental observation in mechanical platforms to the best of our knowledge.

We now further investigate the overtaking behavior by conducting numerical simulations for a longer chain composed of 50 unit cells (see Fig. 5A and its inset for the magnified view). Upon the compressive impact, the TCO chain initially experiences a compression, followed by the tension resulting from the reaction of the TCO cells. In this case, the speed of the initial (i.e., leading) compressive wave is bounded by the system’s sound speed due to its effectively linear nature. However, the following rarefaction pulse is predicted to be supersonic by the KdV analysis. Thus, we would expect the rarefaction pulse to overtake the compressive linear waves (see the “Overtaking behavior” section in the Supplementary Materials for details). We observe this overtaking phenomenon around the fifth cell location (see inset of Fig. 5A), thereby verifying the higher speed of the rarefaction solitary wave than that of the compressive wave pattern.

Fig. 5 Wave speed analysis.

(A) Space-time contour plot of the strain wave for the numerical simulation conducted on the longer chain composed of 50 TCO unit cells. Magnified view of the overtaking moment is shown in the right inset. (B) Trajectory of the rarefaction solitary wave (denoted by the blue markers) and the maximum compressive strain wave (red markers) shows the overtaking behavior of the rarefaction solitary wave. The green line indicates the analytical prediction from the KdV equation. (C) Wave speed of the rarefaction solitary wave is higher than the speed of sound of the medium, which means supersonic behavior.

To quantify the speeds, we extract the time and location of the maximum peak of the rarefaction solitary and compressive strain waves, as shown in Fig. 5B. As marked by the circle in the figure, we detect the overtaking moment around t = 0.11 s when the trajectories of solitary wave and maximum compression peak are crossing. The slope of each curve represents the wave speed; therefore, the wave speed of the rarefaction solitary waves can be calculated from the experimental and numerical results. Figure 5C shows the wave speed calculated from the experiment (20 units) and simulation (50 units), and the result for the rarefaction solitary wave is compared with the prediction from the KdV solution (see the “Continuum limit” section in the Supplementary Materials for the details about the KdV and sound speed calculations). We observe that the experimental and numerical results corroborate the analytical prediction, confirming the supersonic nature of the rarefaction pulse (see the “Overtaking behavior” section in the Supplementary Materials for more details about the overtaking mechanisms in various scenarios).


We have studied experimentally, numerically, and analytically a remarkable example of nonlinear wave propagation in mechanical metamaterials made of volumetric origami cells. We found that the TCO-based metamaterials exhibit the rarefaction solitary wave, which features tensile strains and propagates ahead of the initial compressive strain despite the application of external compressive impact. The KdV-based theoretical analysis provided us with an excellent qualitative handle for understanding the relevant dynamics, and the numerical computations corroborated the experimental observations. Also, the initial compressive strain is attenuated significantly, which can be highly beneficial for impact mitigation applications. While we focused on the monoatomic TCO unit cells with strain-softening behavior in this study, the origami-based system has great potential for supporting rich wave dynamics by introducing heterogeneous elements (e.g., hardening, multistable, and impurity components) in the chain. Also, the findings in this 1D setting can be further extended to multidimensions in a modular way. We believe that this architecture of volumetric origami cells can be used as a versatile building block for a wide range of applications such as impact/shock mitigation, vibration filtering, and energy harvesting.


Prototype fabrication

We used construction paper sheets (Strathmore 500 Series 3-PLY BRISTOL; thickness of the paper is 0.5 mm) for the origami part (see Fig. 1B) and extruded acrylic sheets (thickness of 1.6 mm) for the interfacial polygon (insets of Fig. 1, C and D). These two materials were tailored by a laser cutting machine (VLS 4.6, Universal Laser Systems), as shown in fig. S1A. For the crease of the TCO unit cell, we designed customized crease lines based on the compliant mechanisms (fig. S1B). By using this crease pattern, each crease line shows enhanced fatigue-resistant property, providing repeatable and consistent folding behavior. It should also be noted that this crease pattern (e.g., length, width, and distance of cut slits) affects the compressive force-displacement behavior of the TCO cell significantly. The weight and moment of inertia of the interfacial polygon including a sleeve bearing and markers are 49.4 g and 5.32 ×10−5 kg m2, respectively. The weight of the triangulated facets made of paper is only 3.6 g. Therefore, the inertia of the side facets is negligible compared to the polygon. The assembling process of a single TCO cell is shown in movie S1.

Compression test on the TCO unit cell

The static compression was conducted in a customized testing environment. The unit cell prototype was placed vertically on the load frame (fig. S2), where the bottom part of the unit cell was fixed. The top end of the unit cell was connected to the load cell (LUX-B-50 N-ID, Kyowa) through a sleeve bearing made of polytetrafluoroethylene (PTFE) so that the top polygon was constrained in the axial direction but could rotate freely around the longitudinal axis. The top part holding the load cell was actuated by the motor-driven linear stage (BiSlider, Velmex). The displacement was measured by a noncontact laser sensor (CMOS Multi-Function Analog Laser Sensor IL-065, Keyence).

Preconditioning process on the TCO unit cell

Each TCO prototype was assembled by hand; therefore, the folding behavior, specifically force-displacement relationship, varies among the samples. Similarly, the repeatability of folding/unfolding behavior of each unit cell is also a critical factor to ensure reliable performance of the whole origami system. To address this issue, we conducted fatigue tests on all unit cells by using the load frame explained above. We applied 200 cycles of controlled displacement from −3 mm (tension) to 15 mm (compression) at 6 mm/s, and we measured the corresponding force-displacement relationship. Figure S3A shows the loading and unloading curves from the first cycle (i.e., newly fabricated state without experiencing axial force) to the 200th cycle. One of the noticeable features is that the first cycle shows the significant maximum force peak during loading (denoted by the black curve in fig. S3A). Then, as the number of cycles increases, each curve traces a similar path.

To analyze this fatigue behavior, we considered the area bounded by the loading and unloading curves for each hysteresis loop, and we plotted the area as a function of the number of cycles as shown in fig. S3B. The area enclosed by the hysteresis loop is associated with energy dissipation. For the initial 50 cycles, the energy dissipation decreases significantly, and then the slope of the energy dissipation curve acquires a relatively small value. Hence, to ensure homogeneous and repeatable characteristics of the TCO chain for the dynamic test, we applied this cyclic preconditioning process (200 cycles, displacement of −3 to 15 mm) to all unit cells and selected those with highly consistent and uniform folding behavior.

Dynamic test on the chain of the TCO unit cells

We conducted the dynamic test on the chain of the TCO unit cells by applying compressive impact to the system. For the chain of the TCO unit cells (Fig. 1C), neighboring unit cells were connected mechanically by using M3 stainless steel screws and hex nuts. Also, a flanged sleeve bearing made of PTFE was embedded in the center of the interfacial polygon, through which the stainless steel shaft (diameter of 4.76 mm) was inserted to align the unit cell. In this way, a chain of TCO cells will have free axial and rotational motions while restricting the bending of the chain.

The left end of the chain (first unit cell) was connected to the shaker (LDS V406 M4-CE, Brüel & Kjær) through the customized attachment (see the inset of Fig. 3A) with a sleeve bearing (PTFE), which attaches the leftmost polygon of the cell to the shaker but allows free rotational motions. The right end of the chain (20th unit cell) was fixed to the rigid wall. The shaker was excited by a single-step voltage to apply compressive impact to the system.

To capture the dynamic folding/unfolding motion of the chain, we developed a customized noncontact DIC technique by using Python. To track the motion of each interfacial polygon, we used the color and shape as the features of a target marker. We attached a spherical marker to each corner of the polygons. For the color of the marker, fluorescent green was used to distinguish the marker of the polygons from the other objects. Figure S4A is the original image from the GoPro camera. With the mask based on the fluorescent green, we extracted the marker color area from the raw digital image (see the lower inset of fig. S4A). Last, we identified the marker from this filtered image by examining the shape of the extracted area and determined the 3D coordinate of each marker based on the triangulation method. By using three pairs of the action cameras, we split the field of view horizontally by three and captured the axial and rotational motions of the polygons along the longitudinal axis based on the stereo vision. The six cameras were calibrated before the measurements, and our in-house code processed the rectified images to obtain the 3D coordinate information. Figure S4B shows the axial displacement change of the first interfacial polygon mounted on the shaker attachment. For numerical simulations, we fed this experimentally obtained displacement data into the equation of the first TCO element’s motion.


Supplementary material for this article is available at

Supplementary Text

Fig. S1. Fabrication of the TCO unit cell.

Fig. S2. Compression test on the TCO unit cell.

Fig. S3. Fatigue property of the TCO single unit cell.

Fig. S4. DIC technique to measure the axial displacement and the rotational angle of each polygon.

Fig. S5. Surface plot of the elastic potential energy.

Fig. S6. Effect of damping factor ν.

Fig. S7. Comparison between the numerical simulations on the TCO chain with/without damping effect.

Fig. S8. Comparison between the numerical simulations on the nonlinear/linearized TCO chain.

Fig. S9. Numerical simulation with the application of the initial impact velocity to the first unit.

Table S1. Numerical constants used in the numerical simulation and analytics.

Movie S1. Fabrication of the TCO unit cell.

Movie S2. Experimental demonstration of the rarefaction solitary wave.

Movie S3. 3D reconstruction of the TCO chain from the experimental result.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank S.-g. Ahn and F. O’Brien of the School of Art + Art History + Design at the University of Washington for technical support. J.Y. and P.G.K. also acknowledge the support from the AFOSR (FA9550-17-1-0114). Funding: This material is based on work supported by the NSF under grant no. CAREER-1553202 (J.Y.), DMS-1615037 (C.C.), and DMS-1809074 (P.G.K.). J.Y. and H.Y. are grateful for the support from the ONR (N000141410388) and the Washington Research Foundation. Author contributions: H.Y. and J.Y. conceived the research topic. H.Y. and Y.M. conducted the design, fabrication, and testing of the prototypes. H.Y. conducted the numerical simulations. H.Y., E.G.C., C.C., and P.G.K. created the analytical model and interpreted the computational results. J.Y., C.C., and P.G.K. provided guidance throughout the research. All authors contributed to writing and editing the manuscript. Competing interests: H.Y. and J.Y. are inventors on a patent application related to this work filed by the University of Washington (patent application no. 62/562,596, filed on 25 September 2017). The other authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article