Mechanical cloak via data-driven aperiodic metamaterial design

As shown in Fig. 1A, we focus on the design of a mechanical cloak, ${\mathrm{\Omega }}_c,$ around voids within the structure, that allows it to retain the displacement field in the surrounding region, ${\mathrm{\Omega }}_s,$ as if there were no voids. As a result, the voids are “mechanically invisible” for observers located at any point within ${\mathrm{\Omega }}_s,$ including both the boundaries and the insides. While our method can also accommodate cases that focus only on the boundaries, we choose to adopt this more general and challenging setting, as in most existing studies (24, 26–28), to demonstrate the flexibility of the proposed method. Practically, with this elastostatic cloak, voids could be intentionally introduced within a structure to realize special functions, e.g., to conceal underground tunnels, holes for wires, and cavities to hide objects, without affecting the structural integrity and functionality of the original structure. For instance, the cloak can be used in soft robots to keep the designed distorted postures unaffected by the voids and tunnels to place sensors or driving cables. To achieve this, as shown in Fig. 1B, we consider metamaterials within the structure as homogenized continua and then concurrently optimize the topology of the cloak, ${\mathrm{\Omega }}_c,$ as well as the spatial distribution of the effective properties. Mathematically, we optimize the independent entries of the stiffness tensor, within the cloak ${\mathrm{\Omega }}_c,$ to minimize the void-induced distortion of the displacement field in the surrounding region ${\mathrm{\Omega }}_s$ (see SI Appendix for more details). For simplicity, we consider only unit cells with orthotropic symmetry in this study, whose stiffness tensor is calculated through energy-based homogenization and has four independent entries $C_{11}, C_{12}, C_{22},$ and $C_{33}$ for two-dimensional (2D) design (in Voigt notation) (32, 33). To ensure the fabricability of the optimized materials, the properties are constrained within the properties space of a large and dense precomputed unit-cell database, during the optimization process (Fig. 1D and SI Appendix, Fig. S1). This precomputed database is generated by combining unit-cell topology optimization and a sequential stochastic shape perturbations algorithm to achieve shape and property diversities (see SI Appendix for more details and access to the dataset) (34, 35). Unit cells with isolated pixels, checkerboard patterns, or features smaller than the predefined minimum length scale (0.5 mm) are filtered out to ensure the manufacturing feasibility. As demonstrated in SI Appendix, Fig. S2, the change of the mechanical properties in this precomputed database, i.e., stiffness tensor, is induced by the change of the unit-cell shapes only, without any modification to the constituent materials.

After the optimization, the cloak within the cloaking boundary is filled by an aperiodic tessellation of unit cells, selected to achieve the optimized properties distribution (Fig. 1C shows a printed assembled structure). Specifically, for each location in the cloak, we calculate the infinite norm of the difference between the optimized properties ($C_{11}, C_{12}, C_{22},$ and $C_{33}$) and properties of each unit cell in the precomputed database, which is denoted as ${\theta }_i.$ A set of unit cells is then selected from the database as candidates with the smallest ${\theta }_i$ values, whose properties are closest to the optimized local material properties. A tiling optimization is then performed to select the optimal unit cell from the candidate set in each location, to ensure good geometrical and mechanical compatibility between adjacent unit cells, as shown in Fig. 1E. In the tiling optimization, we use the ratio of the connected region versus the whole shared boundary to measure geometrical incompatibility. To reduce the possible deviation in the stiffness of a unit cell caused by the presence of a neighboring unit cell, we introduce the relative sum of stress difference on the shared boundary under the unit strain field to monitor the mechanical incompatibility. As a result, assembling a well-connected cloak with designed properties distribution can be achieved by selecting the best unit cell from the candidate set of each location to minimize the sum of properties difference ${\theta }_i,$ geometrical incompatibility measures, and mechanical incompatibility measures over the whole structure. This tiling optimization is formulated as an energy-minimization problem on a grid-like graph and solved efficiently by a dual decomposition method (see SI Appendix for more details) (36, 37). After the tiling optimization, the assembled structure with a cloak is obtained (Fig. 1C).

To demonstrate our approach, we first focus on the design of a cloak for a circular void, fixing the shape and size of the cloak. As shown in Fig. 2A, the reference structure (without voids and cloak) is composed of $30\times 30$ periodically tessellated base cells, each of size $5\hspace{0.17em}\times 5\hspace{0.17em}\text{mm},$ whose constituent material has Young’s modulus $E=1.20\hspace{0.17em}\text{GPa}$ and Poisson’s ratio $\nu =\mathrm{0.35.}$ To facilitate displacement tracking in experimental measurement, a simple four-rod cubic lattice is chosen as the base cell (shown in Fig. 2 A, Inset). Its homogenized elastic properties, calculated through energy-based homogenization, are $C_{11}=171.55\hspace{0.17em}\text{MPa}, C_{12}=69.68\hspace{0.17em}\text{MPa},C_{22}=171.55\hspace{0.17em}\text{MPa},$ and $C_{33}=62.07\hspace{0.17em}\text{MPa}.$ Each unit cell is discretized into a $50\times 50$ plane stress Q4 element. Finite-element simulations are then performed to calculate the displacement field of the reference structure under different boundary conditions (see SI Appendix for more details). We first model the response of the structure to a compressive load $u_{bc}=0.5\hspace{0.17em}\text{mm} (0.67\%\hspace{0.17em}\text{average} \text{strain})$ in the horizontal direction, while keeping the top and bottom edges free (we refer to this configuration as a “displacement-free” boundary condition). The response of the material is assumed to remain linear elastic during the loading process. The resultant displacement fields in the homogeneous material, in both the x and y directions, are shown in Fig. 2 B and C, exhibiting a uniform transition from one end to the other. The presence of a circular void at the center of the structure (i.e., in the “voided structure”; Fig. 2F) causes a large distortion of the displacement field (Fig. 2 G and H), with significant shrinkage of the circular void. To control this distortion, we introduce a ring-shaped cloak, ${\mathrm{\Omega }}_c,$ around the void, whose outer radius is twice the radius of the central void, as commonly used in the literature where the cloaking boundary is predefined (24). Later, we also present results when the cloaking boundary is not predefined but instead designed. The homogeneous medium outside the cloak in the structure, filled by periodic base cells as in the reference structure, is denoted as ${\mathrm{\Omega }}_s.$ Using the proposed method, the cloak ${\mathrm{\Omega }}_c$ is filled with aperiodic unit cells selected from the database to minimize the void-induced distortion of the displacement fields within the region ${\mathrm{\Omega }}_s.$ We refer to this structure with a cloak surrounding the void as the “cloaked structure,” which is shown in Fig. 2K. With our tiling optimization, neighboring unit cells retain good connections between each other, even though their properties and geometry vary. It is noted that the proposed optimization method allocates unit cells in a way that gradually increases the stiffness from the outer boundary to the inner one of the cloak ${\mathrm{\Omega }}_c$ (Fig. 2 K and P). This gradient of stiffness compensates for the lack of structural cohesion due to the presence of the void, which is in line with existing mechanical cloak designs. Successful cloaking is evident when comparing the displacement fields of the reference structure, in Fig. 2 B and C, with the one in ${\mathrm{\Omega }}_s$ of the cloaked structure, in Fig. 2 L and M. The large distortion of the central void is mitigated by the optimized cloak, displaying a gradual and uniform transition of displacement field along both x and y directions in the surrounding region ${\mathrm{\Omega }}_s$ as that of the reference structure (i.e., without the void).

As a quantitative measure for the distortion of the displacement field, we opt for the commonly used relative displacement difference $\mathrm{\Delta }$ formulated aswhere ${\overrightarrow{\mathit{u}}}_i$ represents the nodal displacements of the finite elements in the cloaked (Fig. 2 L and M) or voided structure (Fig. 2 G and H), and ${\overrightarrow{\mathit{u}}}_{0,i}$ represents the nodal displacements of the reference structure (Fig. 2 B and C). The summation operator performs the sum over all the nodes within the surrounding region ${\mathrm{\Omega }}_s.$ The lower the relative displacement difference $\mathrm{\Delta },$ the smaller is the overall distortion of the displacement field from the reference state. With the designed cloak, the relative difference calculated from the numerical simulations is reduced from 17.20% in the voided structure to 4.00% in the cloaked structure, demonstrating a good cloaking performance.

To further validate the findings, we perform experiments on three-dimensional (3D)-printed structures, fabricated to reproduce 1) the reference structure, 2) the voided structure, and 3) the cloaked structure (Fig. 2Q). Considering the symmetry of the structure and the camera’s field of view, we measure displacements from the upper left corner of each structure (Fig. 2D). Eq. 1 is used to calculate the relative displacement difference $\mathrm{\Delta }$ in the experiment by inserting the measured displacements of the lattice points (centers of base cells) within the region ${\mathrm{\Omega }}_s.$ Overall, the simulated and experimental displacements match well in terms of both displacement distribution (compare, e.g., Fig. 2 A and B with Fig. 2 D and E or Fig. 2 G and H with Fig. 2 I and J) and the relative displacement difference (reduced from 15.16% in the voided structure to 6.86% in the cloaked structure).

The results shown in Fig. 2 are obtained under the displacement-free boundary condition. Nevertheless, our method is general and valid for various other boundary conditions. We demonstrate this by performing additional optimizations and numerical simulations with pressure-free, pressure-sliding, dilating, and shearing boundary conditions (Fig. 3 A–F and Movie S1). For the pressure-free boundary condition, we apply constant compressive pressure $p_{bc}=95.67\hspace{0.17em}\text{KPa}$ on the left and right edges, to keep the average strain around 1% and ensure that the structure remains in the linear elastic deformation range. The proposed optimization method generates a cloaked structure shown in Fig. 3C, with a similar stiffness gradient as observed in the design obtained for the previous displacement-free boundary condition (Fig. 2 K, P, and Q). As shown in the first two bar groups in Fig. 3B, the relative displacement difference of the voided structure is much higher for the pressure-free boundary condition (102.3%) than that of the previous displacement-free boundary condition (17.2%). This is because, in contrast to the displacement-free boundary condition, the pressure-free boundary condition allows nonuniform displacements on the left and right boundaries, which encourages the shrinkage of the void after the loading (SI Appendix, Fig. S6). Nevertheless, the cloaked structure obtained can still effectively reduce the distortion caused by the void to a $\mathrm{\Delta }=12.8\%.$ Adding an extra sliding boundary condition on the top and bottom in the pressure-free boundary condition leads to the pressure-sliding boundary condition. Under this condition, the constraints imposed by the sliding support compress the structure against Poisson’s-effect–induced expansion in the vertical direction. As a result, compared with the previous single-directional loading cases (displacement-free and pressure-free), the joint effects of the support reactions and imposed pressure enable the void-induced distortion to cover a larger area in the structure (SI Appendix, Fig. S7). Therefore, the optimized cloak contains more unit cells with high stiffness (Fig. 3D and SI Appendix, Fig. S7), to balance the strains and suppress the shrinkage of the void. The expansion of the high stiffness region is even more obvious under the dilating boundary condition (Fig. 3E). This is because, under this boundary condition, the affected area of the distortion in the voided structure is even larger with the uniformly distributed stretching force exerted from all four edges (SI Appendix, Fig. S8). Similarly, for the shearing boundary condition, large void-induced distortions of the displacement field mainly distribute in the right half of the structure (SI Appendix, Fig. S9), corresponding to regions with stiffer unit cells in the cloak (Fig. 3F). Nevertheless, neighboring unit cells in the cloak remain well connected for these three boundary conditions with our tiling optimization, achieving excellent cloaking performance with low relative errors (<9%), as shown in the third through fifth bar groups in Fig. 3B. Note that all the cloak designs shown in this study remain effective for other magnitudes of compressive pressure or displacement imposed on the boundary, as long as the deformation remains in the regime of linear elasticity. This is because, in linear elasticity, the relative displacement difference, Δ, does not depend on the absolute magnitude of the loading.

While most existing design methods for mechanical cloaks confine the base cell in the surrounding region to bimode or simple parameterized lattice, the proposed method can systematically accommodate various free-formed base cells. As an example, we use a different type of unit cell with more complex geometry as the base cell and consider a circular void under the pressure-free boundary condition (Fig. 3G). The homogenized stiffness tensor of this base cell yields $C_{11}=123.02\hspace{0.17em}\text{MPa}, C_{12}=48.02\hspace{0.17em}\text{MPa},C_{22}=123.02\hspace{0.17em}\text{MPa,}$ and $C_{33}=13.50\hspace{0.17em}\text{MPa}.$ As shown by the simulated displacement field (SI Appendix, Fig. S10 and Movie S2) and the large relative displacement difference of the voided structure (Fig. 3B), this type of compliant base cell leads to a significant amount of distortion of the displacement field. Yet, the proposed method can still effectively reduce the relative error from 139.1% to 19.6% with the optimized cloak. Moreover, it should be noted that the proposed method can also design the same cloak to accommodate various loading cases by aggregating individual cloaking performance into an overall objective function. For example, we obtain a cloak that achieves excellent cloaking performance (Δ < 15%) for constant distributed loading imposed from any angle (SI Appendix, Figs. S11 and S12), by using a p-norm–based aggregated objective function in the optimization (see SI Appendix for more details).

Unlike design methods based on the material-parameter transform or direct lattice transform, the proposed method can also be applied to create cloaks around voids with arbitrary shapes (Fig. 3 H–K; SI Appendix, Fig. S13; and Movie S3), due to either functional or aesthetic needs. We demonstrate this capability by designing cloaks for different voids, applying the same pressure-free boundary condition and predefined cloak topology, as in Fig. 2. Similar to our previous designs for a circular void, the stiffness of the unit cells within these optimized cloaks (Fig. 3 H–K) exhibits a decreasing trend starting from the edge of the void to the outer boundary of the cloak. From Fig. 3B and SI Appendix, Fig. S13, it can be noted that the various optimized cloaks reduce distortion of the displacement fields induced by the voids of various shapes, keeping relative displacement difference at a considerably low level (∼6.8% to 13.5%). It is interesting to note that the cloaked structures for inclusions with high aspect ratios and sharp corners (Fig. 3 H–J) can achieve better cloaking performance than that of circular and Mickey-shaped inclusions (Fig. 3 C and K). This might be due to a smaller ratio of inclusion to the cloaking region for structures in Fig. 3 H–J than that for structures in Fig. 3 C and K. It also demonstrates the applicability of the proposed method for inclusions with less regular shapes. It should be pointed out that, while the circular void in Fig. 3C could geometrically enclose voids of various shapes in Fig. 3 H–J, the solution is less optimal or efficient compared to the cloaking design optimized specifically for a desired shape of voids.

In all results discussed so far, a predefined cloaking boundary of the cloak is used and remains unchanged during the design process. This is a requirement similar to most existing methods (24–26, 28). We present next that our proposed method can also simultaneously design the topology, i.e., free-formed geometry that allows topological changes, of the cloaking region ${\mathrm{\Omega }}_c$ and the property distribution within it, adapting to different geometries of voids. To demonstrate the benefits of not predefining the cloaking boundary, based on the Mickey-shaped void, we limit the area of the cloak to be equal to that of the previously used circular one, while concurrently optimizing the topology and its property distribution to minimize the relative displacement difference (Movie S4). The resulting cloaked structure and its displacement field are shown in Fig. 3L and SI Appendix, Fig. S13, respectively. We observe that the cloaking boundary conforms to the contour of the Mickey-shaped void, expanding in the upper half and shrinking in the lower half. As shown in Fig. 3B, while the relative displacement differences are similar for the voided structures with the predefined (88.3%) and optimized topologies of cloaks (86.7%), the cloak with optimized topology can achieve a much smaller relative error (8.1%) than the cloak with predefined topology (13.5%).

To further demonstrate the versatility of the proposed method, we design a cloak for two voids arranged in a butterfly shape, subjected to the displacement-free boundary condition, as shown in Fig. 4. The specialty of this example lies in that there are two neighboring irregular voids (Fig. 4F), which cannot be decomposed into simple ellipses or squares. The existence of multiple voids and the mechanical interaction between them poses an additional challenge to the adaptive design of mechanical cloaks. As shown in Fig. 4 G–J and Movie S5, both the simulated and experimental results suggest that the butterfly-shaped voids distort the displacement field significantly, especially along the y direction. The cloaked structure designed by the proposed method successfully reduced the relative difference Δ nearly by half for both calculated and experimental results, leading to a displacement field (Fig. 4 L–O) close enough to that of the reference structure (Fig. 4 B–E).

In this study, we have developed a data-driven method for the design of mechanical cloaks, powered by topology optimization. This method first optimizes the shape and topology of the cloak and the material property distributions within the cloak and then selects optimal unit cells from a precomputed database of material microstructures to fill the cloak through an efficient optimization-based tiling process. Using multiple examples with various numbers and shapes of voids, boundary conditions, and base cells, we demonstrate that our method achieves excellent performance in mechanical cloaking, verified both numerically and experimentally. Compared to existing approaches that are applicable only to a fixed but restricted topology of cloaks, our method is capable of concurrently designing the topology and property distribution of the cloak, balancing both efficiency and versatility. We observe that the optimized topology of the cloak generally conforms to the contour of different geometries of voids, providing a better cloaking performance compared to cloak designs with fixed topology. While we focus only on voids in this study, the proposed method can be readily applied to the cloak design for solid inclusions with stiffness lower than that of the base cell, which induces smaller loss of structural integrity and thus more easily achieves good cloaking performance. Our method can also handle solid inclusions with a higher stiffness than that of the base cell but will require more compliant unit cells to be included in the database. The data-driven approach is powerful in terms of its capability to swiftly respond to new design scenarios with limited computational cost, by taking advantage of the rich unit-cell database that is readily available. In the preparation of the databases, we can also easily take into account manufacturability and other resource restrictions.

The design of each unit cell in the database is pixelated and represented by a $50\times 50$ binary matrix. We follow the common practice in multiscale topology optimization to use a corresponding $50\times 50$ four-node quadrilateral plane stress finite-element mesh for each unit cell in the full structure. Because the maximal strain of the assembled structure is ∼1% to 2%, the finite-element analysis is carried out in MATLAB and ABAQUS under the linear elasticity and plane stress assumptions. The good agreement between the numerical calculation and the experimental measurement validates the effectiveness of the calculation method.

We fabricated all the structures tested using a Stratasys Connex Objet 500 3D printer. We subjected the structures to compression using an Instron E3000 machine. While the structures were being loaded, we captured a sequence of images using a Nikon D750 camera mounted with a Nikkor 200-mm f/4D IF-ED lens. Digital image correlation (DIC) on the captured images was performed with a global DIC code designed specifically for two-dimensional materials with microstructures (38).

More details on the design methods and experiments can be found in SI Appendix.

We are grateful for support from the NSF Cyberinfrastructure for Sustained Scientific Innovation program (Grant OAC 1835782) and Center for Hierarchical Materials Design (National Institute of Standards and Technology [NIST] 70NANB19H005). L.W. acknowledges support from the Zhiyuan Honors Program for Graduate Students of Shanghai Jiao Tong University for his predoctoral visiting study at Northwestern University.