Flexible mechanical metamaterials are a class of structures with unique geometric features at the microstructural level designed to lead to uncommon properties governed by their nonlinear behaviors [1]. One prominent class of flexible mechanical metamaterials is that of an array of identical pores embedded in an elastomeric sheet [2,3]. Such cellular metamaterials have been shown to be endowed with exotic properties, such as auxeticity [4], programmability [5,6], and tunable wave filtering [7]. These functionalities can be tuned by controlling the exact geometry of the individual pores [8], as well as their arrangement [9]. However, most studies to date have focused on periodic architectures with spatially homogeneous features. On the other hand, aperiodicity has recently emerged as a powerful platform to realize metamaterials with enhanced functionality [10][11][12]. Since the design space is extremely large in aperiodic metamaterials, it is not sufficient to use experiments or direct simulations to discover metamaterials with new properties.
The design of flexible aperiodic metamaterials with desired functionality requires efficient optimization algorithms. While gradient-free algorithms have been successfully used to explore the complex energy landscape of these nonlinear systems [13], they become inefficient as the number of parameters describing the structure increases. On the other hand, gradient-based optimization methods have proven successful in identifying structures with target behavior in the linear regime [14,15], but their ⇤ Corresponding author. E-mail address: bertoldi@seas.harvard.edu (K. Bertoldi). application to the inverse design of nonlinear flexible metamaterials is still nascent [16][17][18]. The success of gradient-based methods to generate flexible metamaterials with target responses ultimately depends on a well-defined objective, the existence of ''good-enough'' nearby local minima (as it is difficult to isolate the objective's global minimum), and the ability to generate new designs that do not violate any imposed constraints.
Two classes of gradient-based methods have been widely used to optimize the response of metamaterials: topology optimization [14] and shape optimization [19]. Topology optimization is a mathematical method that optimizes material layout within a given domain, taking into account specific design constraints. The design variables in topology optimization are continuous fields that specify the presence or absence of material in the domain. Shape optimization, on the other hand, focuses on finding the optimal shape or geometry of a domain with a fixed topology, such as a solid object with embedded a predetermined distribution of holes. In shape optimization, the design variables involve boundary variations rather than the material layout.
Here, we use shape optimization to identify soft, highly deformable, connected, finite, two-dimensional structure with specified stress-strain responses or effective Poisson's ratios. Since topology optimization often fails to provide connected geometries for non-periodic compliant structures (see Appendix), by restricting the space of admissible designs to a fixed topology, we eliminate the possibility of disconnected domains. We start by describing the basic ingredients needed to solve a shape optimization problem. Then, we apply the higher-order moving mesh method to design soft cellular structures with target responses under tension, compression, and shear. We find that all identified solutions are aperiodic, confirming aperiodicity as a promising platform for the design of cellular metamaterials capable of supporting a wide range of mechanical responses in the nonlinear regime.
The ability to withstand large deformations and revert to the undeformed state in a perfectly elastic manner is key to the design of flexible mechanical metamaterials. As a result, we consider an arbitrary domain ⌦ of elastomeric material. The behavior of such an elastic material can be readily described by introducing a strain energy functional that depends solely on the local material deformation described by the displacement field u = (u x , u y ) T . Here, we capture the material response using a nearly incompressible neo-Hookean model with strain energy density function
where F = ru + I is the deformation gradient, µ is the shear modulus, and is the first Lamé coefficient. The neo-Hookean model is a widely used constitutive model that well describes the mechanical behavior of the elastomeric materials typically used to realize flexible mechanical metamaterials [4]. We note that our approach is not restricted to neo-Hookean, and could be applied to a broad class of hyperelastic constitutive models.
Equilibrium of the structure is ensured by finding the displacement field u 2 V that minimizes the total energy ⇧(u), where V is the space of kinematically admissible displacement fields. The total energy is defined as
where t is the traction applied on a portion of the boundary N and b represents a body force. This amounts to solving the variational problem @⇧(u; v) = d⇧ (u + ✏v)
and can be equivalently written as
In this study we use the finite element (FE) method to solve Eq. ( 4) for different loading conditions and domains ⌦ and apply shape optimization to identify the domain shape that leads to target mechanical responses.
Shape optimization seeks to identify the optimal shape of a domain ⌦ that minimizes a shape functional J : U ad ! R, where U ad denotes the set of admissible domains. In the context of this study, J evaluates the mismatch between a desired target structural properties (e.g. target stress-strain curve and target tunable Poisson's ratio) and the response of the structure. Generally, the shape optimization problems considered in this study can be written as find ⌦ ⇤ = arg min
where ⌦ ⇤ is the optimal geometry, @⇧(u, v) is defined in Eq. (3) and c i (⌦) denotes the ith additional constraint imposed on the domain. Given the high sensitivity of flexible mechanical metamaterials to small variations in geometric features [8], we expect shape optimization to be a powerful tool to design a system with desired properties. However, these problems are difficult to solve since the dependence of J on the domain is typically non-convex [19] and non-unique.
In this study, we employ a higher-order moving mesh method to optimize the shape of the considered structures [20]. While other strategies have been proposed such as the level set method [21], phase field method [22], and reduced mapped parametrizations [17], the higher-order moving-mesh method offers three main advantages: (i) it enables high boundary resolutions via an arbitrarily refinable conforming mesh, (ii) it enables the exploration of a large class of shapes beyond polytopes, and (iii) it is inherently compatible with standard finite element simulation packages [20]. The higher-order moving-mesh method represents the control space as a set of invertible geometric transformations T k that can be parameterized by the same underlying finite element mesh used to represent the solution u. From a mechanics perspective, the mapping T k can be thought of as a new intermediate transformation of the initial domain ⌦ 0 to a reference configuration ⌦ k (see Fig. 1).
We simulate the response of flexible mechanical metamaterials under different loading conditions by solving Eq. ( 4) using the open-source FE solver firedrake [23]. In all of our analysis, we use plane strain conditions and discretize the models with quadratic triangular elements unless otherwise specified. We employ adaptive continuation to improve convergence of the forward problem [24]. Adaptive continuation splits up the desired interval into smaller domains and will refine the step size if the previous state does not generate a good-enough starting point for the current iterate. In problems where adaptive continuation is not sufficient, we instead perform a dynamic step using a backward Euler time discretization to generate a guess for the current load increment. We also make use of an efficient checkpointing routine. The checkpointing routine records the load-state pairs of the previous shape iteration and uses that as an initial guess when we solve our nonlinear equations via Newton's method. By checkpointing we avoid recomputing the entire trajectory via adaptive continuation or a dynamic step and reduce the number of forward solves especially in steps where the shape of the domain does not change significantly.
To identify novel geometries leading to target behaviors we use fireshape, a shape optimization library built atop firedrake that is based on the higher-order moving-mesh method [25]. Sensitivities and the update are computed automatically through the fireshape interface. We solve the resulting problem using the augmented Lagrangian method which facilitates the introduction of additional constraints. The unconstrained nonlinear optimization problems that are generated by the augmented Lagrangian are solved using a trust-region algorithm with a quasi-Newton L-BFGS Hessian approximation [26]. We configure the trust-region solver to have an initial trust radius of 10 3 L (where L is the system's characteristic length scale), which prevents large initial update steps that can lead to excessive mesh distortions. An explicit barrier is constructed for problems that fail to converge numerically and for steps that lead to excessive distortion. In cases where a FE simulation fails to converge or the gradient of the shape transformation is below a threshold det(rT k ) < 0.05 (evaluated for each element) we assign a NaN to the objective function and the trust radius. In this case, the trust radius is decreased by a factor of 0.25L and the optimizer computes a new step from the previous iterate. The optimization problem is stopped when either a maximum of 50 iterations or a minimum step size of 10 10 L or a minimum gradient tolerance of 10 6 are reached.
Our starting point is a structure consisting of an elastomeric square block with edges of length L and an embedded 9 ⇥ 9 square array of circular holes (see Fig. 2-a). The holes have radius r and center-to-center distance a = L/10, chosen so that the initial porosity is 0 = ⇡ r 2 /a 2 = 0.6. Further, the structure has the two vertical edges flanked by a column of semicircles and the horizontal ones ending with a strip of solid material of width, a r. As such, the area of the structure covered by elastomeric material is A elas = 0.462L 2 . The behavior of the elastomeric material is captured by a quasi-incompressible neo-Hookean model with strain energy density given by Eq. ( 1), initial shear modulus µ and first Lamé coefficient = 24µ, which results in a Poisson's ratio of ⌫ = 0.48. In Fig. 2-b, -c and-d we report the evolution of the normalized nominal stress as a function of the applied strain for the structure subjected to uniaxial tension, uniaxial compression and shear, respectively. To simulate uniaxial loading conditions we apply a normal displacement on the top and bottom boundaries, u n = 0.5L " n, where " is the applied strain and n denotes the outward unit normal to the edges. During these simulations we monitor (i) the resulting reaction force in vertical direction at the top boundary per unit thickness, F 2 , from which we calculate the nominal stress as S 22 = F 2 /L, and (ii) the displacement in horizontal direction in a region of size L/2 centered at the middle of the structure on the left and right boundaries, L ⌫ and R ⌫ (see region highlighted in blue in Fig. 2-a), from which an effective structural Poisson's ratio is then calculated as
In the shear simulations we apply a tangential displacement on the top and bottom boundaries, u ? = 0.5L"n ? , where n ? is a unit vector obtained by rotating n by 90 in a clockwise direction. We monitor the resulting reaction force in horizontal direction at the top boundary, F 1 , from which we calculate the nominal shear stress as
The results of Fig. 2-b indicate that under uniaxial tension all pores elongate in the direction of the applied strain, leading to a nearly-linear stress-strain curve with S 22 /µ ⇡ 0.94" and a slowly decreasing effective Poisson's ratio ⌫ ⇡ 0.27 0.45". By contrast, when subjected to uniaxial compression the response of the structure is characterized by two distinct regions (Fig. 2c): (i) an initial nearly-linear stress-strain response with stiffness identical to that measured upon uniaxial tension, S 22 /µ ⇡ 0.99", and nearly-constant Poisson's ratio, ⌫ ⇡ 0.28 and (ii) for strains in excess of " cr = 0.025 a stress plateau and a rapidly decreasing effective Poisson's ratio that eventually becomes negative. In the initial linear regime all vertical ligaments separating the circular holes uniformly compress, but the sudden departure from linearity is caused by their buckling, which induces the formation of a pattern of mutually orthogonal elliptical holes [2] (see inset in Fig. 2-c). This pattern of perpendicular elongated holes leads to an effective negative Poisson's ratio, and substantially reduces the overall stiffness of the structure. Finally, for shear we find a nearly-linear stress-strain response with S 12 /µ ⇡ 0.06" (Fig. 2-d).
Because of its rich mechanical behavior the structure presented in Fig. 2 provides a good platform to identify architectures supporting a wide range of mechanical responses. As such, we choose it as the initial domain ⌦ 0 and use shape optimization to find domain shapes that lead to target stress-strain responses and Poisson's ratio evolutions under uniaxial tension, uniaxial compression and shear. In pursuit of a target response, we define the shape functional (objective) J as the squared difference between the computed response f (") and target behavior f targ (x) over the considered range of applied strain " 2 [0, " max ],
where the function f stands for either the nominal stress or the effective Poisson's ratio and the integral is approximated using a composite Simpson's quadrature rule with N = 12 equally spaced control points " i through the quadrature weights w i . Finally, we note that in all our optimization analyses we choose |" max | = 0.1, constrain the top an bottom edges to remain horizontal and of length L and keep the total area of the elastic matrix identical to that of the initial domain (i.e. 0.462L 2 ). When seeking for structures with target stress-strain response, we allow for changes in shape of the boundaries of all the pores. Differently, when looking for target Poisson's ratio evolutions, to facilitate the calculation of the Poisson's ratio, we do not allow for changes in shape of the two columns of semicircular holes that flank the structures.
We aim to identify geometries that under uniaxial tension result in a target stress-strain response
where
denotes the initial stiffness and C targ 2 determines the stiffening behavior. Such geometries are identified by minimizing the shape functional J as defined in Eq. ( 6) with f (") = S 22 (")/µ and f targ (") = S targ 22 (")/µ. Since flexible mechanical metamaterials subjected to large deformations typically exhibit non-linear stress-strain behavior, we first aim at identifying architectures characterized by a purely linear response when uniaxially stretched. More specifically, in Fig. 3 we report results for C targ 2 = 0 and C targ 1 = 0.1, 0.2, 0.3, 0.5, 0.8, 1.0. In Fig. 3-a we show the evolution of the objective function J at each instance it was evaluated, in Fig. 3-b we compare the optimal stress-strain curves to the target ones and in Fig. 3-c we report the optimal geometries in the undeformed configuration (i.e. at " = 0) and at " = 0.1. We find that for all considered targets the optimization algorithm converges after approximately 100 function evaluations and that these solutions closely match their corresponding targets. Further, inspection of the identified optimal geometries reveal that none of them is periodic. As the overall structural stiffness increases, the region of elongated mutually orthogonal holes monotonically shrinks. Further, when the stiffness of the target response approaches or surpass that of the initial domain ⌦ 0 , the optimization algorithm allocates more material to the vertical ligaments, effectively creating thick vertical columns that accommodate the imposed axial deformation and lead to a stiff response.
The geometries identified by the optimization algorithm using the elasticity norm exhibit stress-strain curves that closely match the target ones, but given the extremely large design space we expect many more architectures to exist that result in similar mechanical responses. For example, the response of all designs remains unaltered when they are reflected over the x-axis and the y-axis. However, due to the deterministic nature of the employed algorithm, we are unable to access multiple solutions leading to the same response using the same initial domain and numerical scheme. To identify them we must either change the shape of the initial domain or modify the optimization algorithm. Though, there are no guarantees that either of these methods will identify different structures leading to the same response. In an attempt of identifying a different solution to the optimization problem (while avoiding symmetric solutions) we change the update in the optimization algorithm and choose to calculate it as the Riesz representation of the H 1 norm. In Fig. 4 we focus on structures with C targ 2 = 0.0 and C targ 1 = 0.2 and 0.5 and report results obtained using the H 1 norm. Our results indicate that for C targ 1 = 0.5 the optimization algorithm finds a geometry whose response closely approximates the target. Remarkably, this optimal structure is very different from that reported in Fig. 3-c and it is nearly periodic. All pores have shapes close to the initial circular one and the target stiffness is met by thinning the ligaments in the vertical direction. However, this design strategy fails when considering more compliant targets. For C targ 1 = 0.2 the objective function is hardly reduced since the optimization algorithm attempts to converge to structures with very thin vertical ligaments whose behavior leads to large distortions in the underlying mesh. As described in Section 4, when the mesh is excessively distorted (i.e. for det(rT) < 0.05), we abort the simulation and assign a NaN to the objective function, preventing convergence to a local minimum. ). Similar to the structures presented in Fig. 3, all optimal geometries are non-periodic and comprise a cluster of mutually orthogonal elliptical holes surrounded by nearly-circular pores. However, different from the optimal geometries for linear response, in this case the size of the cluster with elliptical holes is similar for all four geometries and the stiffening rate is controlled by the aspect ratio of the pores. As we move from C targ 2 = 1 to C targ 2 = 5, the aspect ratio of the elliptical pores monotonically decreases. This is because the more elongated the elliptical holes are, the more applied deformation is needed for the ligaments to align along the vertical direction and, therefore, lead to a stiffer response.
To validate the numerical results, we fabricated the optimal structures identified by our optimization algorithm for (C ) = (0.3, 0.0) and (0.2, 3) via a molding approach out of a silicone rubber (Elite Double 32, Zhermack -with µ = 0.262 MPa [27,28]). In Fig. 7-a, we compare the numerical and experimental stress-strain curves for the two structures, whereas in Fig. 7-b we display snapshots taken at " = 0.0 and 0.1.
We find good agreement between simulations and experiments, confirming the validity of our optimization strategy.
Having demonstrated that the optimization algorithm can successfully identify architectures with target stress-strain curves, we next seek to find architectures with a desired effective Poisson's ratio ⌫ over a range of applied strain,
In Fig. 6 we present results for both target constant (with (D ( 0.1, 0.5) the optimization algorithm successfully reduces the objective function J and identifies geometries with effective Poisson's ratio that closely matches the desired one. The optimal structures discovered are again aperiodic and comprise a cluster of mutually orthogonal elliptical pores. Since the geometry of this region resembles that of the rotating squares, which are known to exhibit large negative values of Poisson's ratio [29], the size of the cluster grows as ⌫ is reduced. An exception is represented by the structure with a specified target of ⌫targ = 0.1 0.5". In this case the optimization algorithm identifies another mechanism to achieve large negative values of ⌫. Specifically, it combines a smaller region of elongated ellipses with pores that resemble rotated rounded squares. However, although this geometry leads to an effective Poisson's ratio ⌫ more negative than the other considered in Fig. 6, its response does not match the decreasing and(E 0 1 , E 0 2 , 0.024) (light blue), where E 0 1 = 0.996 and E 0 2 = 0.0167 are the pre-and post-buckling stiffness measured when compressing the initial domain ⌦ 0 . (b) Target response (circular markers) and stress-strain curves for the optimal designs (continuous lines). The black dotted line denote the stress-strain curve of the initial domain ⌦ 0 . (c) Numerical snapshots of the optimal geometries in the undeformed configuration (i.e. at " = 0) and at " = 0.1. (d) Zoom in of the optimal configuration (dark pink) overlayed above the initial domain ⌦ 0 (black).
(For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) target behavior. This discrepancy may be due to either the fact that this specific target behavior is unfeasible or is far from the initial starting configuration.
We now seek to optimize structures with target behavior upon uniaxial compression. This is a more challenging task since instabilities can be triggered upon compression, which can make the underlying system of nonlinear equations poorly conditioned and result in a complex energy landscape, characterized by multiple bifurcations from which stable and unstable solution branches emanate [5].
We aim at identifying structures with target initial stiffness and an approximate critical buckling strain. To this end, we specify a target response of the form
During optimization we do not explicitly prescribe the critical buckling strain " targ cr , since this would require either the introduction of an additional eigenvalue constraint [30] or a formulation constrained by an augmented system of equations [31]. Instead, as described in Section 4, we still use 12 equally spaced points to describe the target and calculate J (see Eq. ( 6)). In Fig. 8 we report results for (E and(E 0 1 , E 0 2 , 0.024) where E 0 1 = 0.996 and E 0 2 = 0.0167 are the pre-and post-buckling stiffness measured when compressing the initial configuration ⌦ 0 . For all three cases the optimization algorithm reduces the objective function J over the iterations (Fig. 8-a) and identifies architectures with responses that closely matches the target ones (Fig. 8-b). All optimal structures reported in Fig. 8-c at first glance appear to be identical to the initial domain, with a square array of circular pores. However, closer inspection of the holes reveals the target behavior is achieved by introducing small fluctuations that move the boundaries of the pores away from a perfect circular shape (Fig. 8-d). Since these underlying aperiodic fluctuations do not suppress the instability, it suggests that the buckling-induced pattern transformation observed in this class of mechanical metamaterials is a robust phenomenon.
Additionally, we can use the optimization algorithm to identify architectures for which instabilities are not triggered upon application of a compressive strain " = 0.1. Towards this end, we again use Eq. (7) to specify a linear stress-strain target response defined by C targ 2 = 0.0 and C targ 1 = 0.2 and 0.4. To our surprise, the optimization algorithm did not converge and could not identify geometries that support the desired responses. In order to improve convergence, we then switched to linear elements and refined the original mesh. While with these changes the algorithm can perform more iterations, it cannot still identify geometries with responses that closely match the target ones (Fig. 9-b). More specifically, for the case of C targ 1 = 0.4 the identified structure nearly converges to the target linear response, whereas for C targ 1 = 0.2 no architecture is found that captures the target response.
Next, we examine the stress-strain curves of the intermediate domains ⌦ k explored by the optimization algorithm when searching for linear behaviors and find that many of them can be captured by Eq. ( 7) with C for these two targets the objective J is greatly reduced over the iterations and that the identified optimal geometries closely match the target responses. We note that both structures are once more aperiodic and comprise domains of mutually orthogonal elliptical pores. In these regions the ligaments are not straight in the undeformed configuration and therefore mainly deform via bending, suppressing the instability. As for the exact mechanisms governing the evolution of the overall structural stiffness, these are hard to deduce. However, the ability of the identified designs to capture a wide range of mechanical behaviors suggests that aperiodic metamaterials may support a richer set of functionalities than their periodic counterparts.
Finally, we focus on tailoring the effective Poisson's ratio under compression. While the effective Poisson's ratio for the initial domain exhibits a sharp transition induced by buckling (Fig. 2-c), here we aim at identifying architectures with constant ⌫ over the considered range of applied compressive strain. More specifically, in Fig. 10 we report results for five constant target Poisson's ratios specified by Eq. ( 8) with D targ 1 = 0 and D targ 0 = 0.2, 0.1, 0.0, 0.1 and 0.2. The snapshots of the optimal structures reported in Fig. 10-c indicate that to achieve the target behaviors the optimization algorithms identifies a mechanism different from that exploited to reach a target Poisson's ratio under tension (see Fig. 6). More specifically, the algorithm converges to structures with several rows of pores next to the horizontal boundaries that transform into a pattern of mutually orthogonal elongated ellipses upon compression. This transformation leads to a significant lateral contraction, but it does not affect the effective Poisson's ratio, since it occurs outside of the region used to evaluate ⌫. Further, it absorbs most of the applied compressive strain, so that the interior region does not have to deform significantly to satisfy the boundary conditions. Finally, since the chosen targets have a lower Poisson's ratio compared to the initial domain in the undeformed configuration, the pores of the interior region are mutually orthogonal ellipses, which aspect ratio that monotonically increases as ⌫ decreases.
We now search for geometries that under shear result in a target linear response
In Fig. 11 we focus on C targ 1 = 0.025, 0.05, 0.08, 0.1. For all prescribed targets the optimization algorithm finds aperiodic and non-intuitive pore arrangements leading to stress-strain responses that closely match the target ones. To achieve the = 0.025 (dark green), 0.05 (teal), 0.08 (light blue), 0.1 (light pink). (b) Target response (circular markers) and stress-strain curves for the optimal designs (continuous lines). The black dotted line denote the stress-strain curve of the initial domain ⌦ 0 . (c) Numerical snapshots of the optimal geometries in the undeformed configuration (i.e. at " = 0) and at " = 0.1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) softest target behavior the numerical algorithm tilts the vertical ligaments in the central part of the structure in direction opposite of that of the applied strain. As a result, these ligaments mainly rotate upon shearing, leading to a soft behavior. By contrast, the vertical ligaments in the stiffest structure are tilted in the direction of the applied strain and generate a diagonal band-like structure. It follows that the applied deformation further stretch them, resulting in a stiff response.
We have used shape optimization to identify porous architectures with target mechanical responses. In contrast to previous studies which have focused on the optimization of periodic unit cells [16,17], here we considered a finite size structures. Remarkably, all optimal solutions identified by our numerical algorithm are aperiodic and non-intuitive. As such, our study suggests that non-periodicity may open new avenues for the design of structures capable of supporting a wide range of mechanical responses.
In this study we focused on an elastomeric metamaterial comprising a square array of holes and subjected to uniaxial tension, uniaxial compression, and shear, but the proposed numerical approach can be extended to flexible structures with different geometries and made of different materials. As an example, in Fig. 12 we consider an elastomeric square block with edges of length L and an embedded triangular array of holes. We choose the holes to be circular in the initial domain ⌦ 0 , with radius r = 0.0442L and center-to-center spacing of L/10, resulting in an area A elas = 0.325L 2 covered by elastomeric material. When subjected to uniaxial tension, this structure has a nearly-linear stress-strain curve with S 22 /µ ⇡ 0.325" (dotted line in Fig. 12-b).
We then seek a domain shape leading to a stiffer linear response, S 22 /µ = 0.4". As shown in Fig. 12-c, the optimization algorithm successfully identifies an aperioidic structure whose stress-strain response closely approximates the target.
Our results indicate that the proposed optimization strategy is capable of identifying optimal geometries for different loading conditions and initial geometries, but for certain targets it may fail to find an appropriate structure. Such failures arise either because these targets are unfeasible, or because the designs separating the initial domain from the local minimum exhibit very small geometric features that cause excessive distortion in the underlying mesh. To overcome this issue, in future endeavours the mesh-distortion constraint can be relaxed in place of a method that adaptively remeshes to preserve an underlying high quality mesh.
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.
In Fig. A.13, we plot the evolution of the strain energy and the overall volume of the structure for both uniaxial tension (left) and shear (right). We find that for both loading conditions the algorithm generates disconnected domains. To the best of the authors' knowledge, the ability to generate connected finite sized aperiodic still remains an open problem in topology optimization. As such, in this study we chose to use shape optimization to guide the design of nonlinear flexible mechanical metamaterials.