Optimization-based design of an elastostatic cloaking device

We present a new method for the design of devices to manipulate the displacement field in Elastic materials. It consists of solving a nonlinear optimization problem where the objective function defines the error in matching a desired displacement field, and the design variables determine the required material distribution within the device. In order to facilitate fabrication, the material at a given point of the device is chosen from a set of predefined materials, giving raise to a discrete optimization problem that is converted into a continuous one using the Discrete Material Optimization technique. The candidate materials maybe simple, isotropic materials, but the device made of them behaves as a whole as a metamaterial, enabling the manipulation of the displacement field in ways that are inconceivable in nature. As an example of application, a device for elastostatic cloaking or unfeelability is designed.

We present a new method for the design of devices to manipulate the displacement field in Elastic materials. It consists of solving a nonlinear optimization problem where the objective function defines the error in matching a desired displacement field, and the design variables determine the required material distribution within the device. In order to facilitate fabrication, the material at a given point of the device is chosen from a set of predefined materials, giving raise to a discrete optimization problem that is converted into a continuous one using the Discrete Material Optimization technique. The candidate materials maybe simple, isotropic materials, but the device made of them behaves as a whole as a metamaterial, enabling the manipulation of the displacement field in ways that are inconceivable in nature. As an example of application, a device for elastostatic cloaking or unfeelability is designed.
Cloaking, that is hiding objects to certain fields, was first considered by Hashin and Shtrikman 1 , who found that spheres with an appropriate coating do not disturb the magnetic flow in the surrounding material.
In 2006, simultaneously and independently, Leonhardt 2 and Pendry et al. 3 introduced the use of conformal mapping to determine the inhomogeneous and anisotropic refractive index for electromagnetic cloaking, giving birth to the approach based on transformation optics for metamaterial design. The inhomogeneity and the anisotropy of the cloaking medium are the result of applying the conformal mapping to a homogeneous medium under uniform electromagnetic field. The required materials may have extreme properties (for instance, negative refraction index). These properties are usually unconceivable in nature, and a composite has to be designed to mimic them; this is the so-called metamaterial.
Applied to Elasticity, Milton et al. 11 showed that transformation optics produces anisotropic Cosserat materials 12 . Later, Norris and Shuvalov 13 demonstrated that the transformation method produces Willis materials 14 , having Cosserat materials as simplified particular cases. Considering the little knowledge on the realization of specific Willis or Cosserat materials, Bückman et al. 15 proposed the "direct lattice" transformation approach for the design of realizable mechanical metamaterials from lattices.
As alternative for the design of realizable metamaterials, we proposed the optimization-based approach for metamaterial design 16 . It consists of solving an optimization problem where the objective function measures the accomplishment of the task (e.g., cloaking) assigned to a metamaterial device and the design variables define the distribution of parameters describing the microstructure of the metamaterial in the device. Like the direct-lattice transformation approach 15 , the optimization-based approach directly prescribes how to fabricate the metamaterial at a point, but it excels the former in allowing not only lattices but any quantitatively characterized material.
However, beyond the difficulty yet impossibility of fabricating a specific anisotropic metamaterial, all the above mentioned approaches produce inhomogeneous metamaterials, which is another major deterrent for the design of practical metamaterial devices.
In previous works 17, 18 , we introduced the idea of considering a device for heat flux manipulation to behave as a metamaterial as a whole but being made of a few simple even isotropic materials. This rebuts the mainstream belief that it is imperative to use anisotropic inhomogeneous metamaterials for the manipulation of macroscopic fields. Now, this way of designing what are maybe the easiest to fabricate devices for the manipulation of macroscopic fields is extended to elastostatic problems. Particularly, its capability is proven by means of the design of a device for elastostatic cloaking or unfeelability using only two isotropic elastic materials.

Methodology
Let Ω be originally a body made of a given arbitrary material, and u 0 be the displacement field in Ω under given tractions t and displacements u in the portions ∂Ω t and ∂Ω u , respectively, of the boundary ∂Ω of Ω. Then, let us assume that an inclusion of a second material, with either much less or more stiffness than the original one, is embedded in the region Ω incl ⊂ Ω, as shown in Fig. 1. The presence of the inclusion sensibly affects the displacement field in Ω, which is now referred to as u.
Using the finite element method (FEM), the displacement u at any point x ∈ Ω is approximated as follows: n N n n 1 nod where ϕ n is the shape function associated to the node n of the finite element mesh of Ω, n = 1, 2, …, N nod , and u n is the (unknown) displacement at this node; for  ∈ u N dim (N dim = 2 for plane strain and plane stress states, d of is the matrix grouping the shape functions, N dof = N dim N nod , and  ∈ U N dof is the vector of nodal displacements, whose components u n = u(x n ) are prescribed for all the nodes x n ∈ ∂Ω u . On the other hand, the unknown components of U are determined as solution of the equilibrium equations: where K and F are the global stiffness matrix and the nodal load vector, respectively, given by T t with B as the strain/displacement matrix and C as the effective elastic moduli. The linear algegraic system of equations (2) is the FEM version of the equilibrium equations for linear elastic solids, whose solution is widely detailed in the literature (see for instance the book of Zienkiewicz and Taylor on the basics of FEM 19 ). Now, to cloak the inclusion Ω incl requires to have u = u 0 at the points located in a certain region Ω cloak where the displacement is to be sensed. To this end, we must design the material inside a certain region Ω dev (from now on, referred to as the "device") surrounding Ω incl . In general, such material by design or metamaterial has a variable microstructure throughout the device Ω dev .
Given Ω divided in finite elements, let the microstructure at each finite element Ω (e) ∈ Ω dev be quantitatively characterized by N par scalar parameters grouped into the vector p (e) ; this means that any effective material property at Ω (e) can be expressed as a function of p (e) , as it is the case of the elastic moduli C = C(p (e) ). The microstructure throughout the device is characterized by the vector P made of all the parameters p (e) of all the N dev finite elements in the device Ω dev , so  ∈ P N var with N var = N par N dev . Consequently, being K = K(P) in the equilibrium equations (2), it immediately follows that u at any x ∈ Ω depends on P, i.e. u = u(x, P).
Then, the cloaking design problem consists of finding P (i.e., the microstructure distribution in Ω dev ) such that u(x, P) = u 0 (x) for all x ∈ Ω cloak . To make such problem amenable to be numerically solved, instead of checking the accomplishment of the cloaking task at all the points x ∈ Ω cloak , let us do it at N check predefined checking points X i ( ) ∈ Ω cloak . So, the discrete form of the cloaking design problem can be stated as: to find  ∈ P N var such that par defines an admissible microstructure, which constrains the search of P to a feasible design set  ⊂ N var  . In general, it will not be possible to exactly accomplish the cloaking task (5) by looking for P in the set . Therefore, in order to obtain an optimal design, we propose to solve the following nonlinear constrained optimization problem: obj where P plays the role of decision or design variables, N var is the number of design variables, and f obj is the objective function, defined as which is the root mean square error (RMSE) in the accomplishment of the cloaking task (5). Using the current approach, i.e. defining the cloaking task as the optimization problem (6), it is still possible to obtain a design for perfect cloaking if the feasible design set  is rich enough. If not, we will obtain a design for which the error in the accomplishment of the task reaches a minimum.
Parametrization of the microstructure. In a broad range of materials, the so-called "quantitatively characterized" materials 20 , the effective material properties depend on a few parameters: the thickness and the orientation of layers in laminates 16 , the density and the irregularity factors in materials with isolated inhomogeneities 20 , the size of the prismatic inclusions in an elastic composite 21 , the fiber orientation in fiber-reinforced polymers 22 , the size of particles or beads in coating of dental implants 23 , etc. These parameters are the components of the vector p (e) at a finite element.
The metamaterials designed using the direct lattice transformation approach 15 are quantitatively characterized materials where p (e) defines the topology of the unit cell of the lattice at the element Ω (e) . Actually, Cosserat materials like those designed using the transformation optics approach can be represented by assuming the components of p (e) to be the effective material properties (or their components if they are tensors), following the free material optimization (FMO) approach 24 . By this way, the current optimization-based approach for metamaterial design is capable of embodying the previous design approaches as particular cases.
After designing a device for thermal cloaking using variable metamaterial distribution 16 , and realizing the difficulty of fabricating it, we proposed 18 that the material at any point of the device should be chosen from a list of previously defined candidate materials. This approach, leading to easy-to-make solutions, is adopted in this work.
Then, after adopting a set of candidate materials consisting of N cand linear elastic solids with known elastic moduli C m , m = 1, 2,…, N cand , the effective elastic moduli C (e) at a finite element Ω (e) ∈ Ω dev is defined by the mixture law:  optimization problem, which is too expensive to solve when the number of design variables is large (as it is the case when there are N cand variables per finite element in a fine enough mesh).
Then, we transform this integer optimization problem into a continuous one, which can be efficiently solved using gradient-based algorithms. To this end, we use the "discrete material optimization" (DMO) technique 25

Results
Let us apply the current methodology for the design of a mechanical cloaking device similar to that one designed by Bückmann et al. 15 using the lattice transformation approach. Given a holed plate Ω made of nylon, compressed under the load 100 kN/m applied at two opposite faces (see Fig. 2), the problem consists of cloaking the hole Ω incl using a coating ring that occupies the region Ω dev .
As candidate materials for building the cloak, we chose two isotropic materials: aluminum, with Young modulus E = 69 GPa and Poisson ration ν = 0.32, and polytetrafluoroethylene (PTFE), with E = 0.5 GPa and ν = 0.4. The stiffness of these materials are sensibly different than that of the base material (nylon has E = 3 GPa and ν = 0.4).
The plate Ω is assumed to be under plane stress conditions and modeled using a mesh of 200 × 200 bilinear square finite elements.
Considering the nylon plate without the hole, the FEM solution for the displacement field, that is u 0 , is that given in Fig. 3(b,c).
For the nylon plate with the hole Ω incl , let us take the previous mesh and then discard the finite elements whose centers lie in Ω incl . In this case, the FEM solution for the displacement is that shown in Fig. 3(e,f), and the error in the accomplishment of the cloaking task is = .
= .  u RMSE 1 845 mm 0 685 max nocloak 0 , with = . u max 2 693 mm 0 . Now, let us design a device Ω dev around the hole Ω incl such that the displacements approach u 0 at all the N check = 19972 nodes outside Ω dev . Using the current approach, we seek to accomplish the cloaking task with a minimum error by solving the optimization problem (6) with the design variables P defining the material at all the N dev = 15084 finite elements whose centers lie in Ω dev . So, the total number of design variables is N var = 2N dev = 30168. In this work, this large nonlinear constrained optimization problem (6) is solved using the interior-point algorithm known as IPOPT 27 . As suggested by Sigmund 28 , a smoothed Heaviside density filter was applied to reduce both checker-board type instabilities and "grey zones" (those where the material is none of the candidates but a composite of them). Figure 3(g) shows the so-computed optimal material distribution in the cloaking device, giving the displacement field depicted in Fig. 3(h,i), with RMSE = 0.00417RMSE nocloak .
To eliminate the grey zones that are observed in the optimal design shown in Fig. 3(g), let us apply the simple a posteriori black-and-white filter proposed by Fachinotti et al. 17 : if the aluminum fraction at an element is greater than a certain threshold φ * ∈ (0, 1), the element is assumed to be fully made of aluminum; otherwise, it is made of PTFE. By this way, the device becomes easier to fabricate at the expense of deteriorating the accomplishment of the cloaking task since the material distribution is not longer optimal. A priori, the natural choice for φ * should be φ * = 0.5, that is equivalent to assume that the material at an element is that having the highest fraction; however, the corresponding device (Fig. 3(j)) poorly performs the cloaking task, with RMSE = 0.5495RMSE nocloak . Actually, the lowest error in the cloaking task is achieved for the quite low threshold φ * = 0.015, giving raise to the device shown in Fig. 3(m), for which RMSE = 0.0610RMSE nocloak . For easier fabrication, let us eliminate the narrow aluminum areas at the top and at the bottom of the device shown in Fig. 3(m), obtaining that shown in Fig. 3(p), for which RMSE = 0.0896RMSE nocloak . This device has an aluminum core with a variable thick coat of PTFE, the whole acting as a compliant mechanism to cloak the hole.
Considering the measure of the error in the task accomplishment as defined by Bückmann et al. 15 , that is we obtain Δ = 0.40% for the optimal device (that of Fig. 3(g)) and Δ = 8.56% for the simplest device (that of Fig. 3(p)). This is a very good cloaking performance compared to that obtained by Bückmann et al. 15 , for whom Δ ≈ 20%. Extension to multiple loads. This methodology can be easily applied to accomplish another tasks by just changing the objective function. Further, to change the involved geometries (either that of the device, the inclusion or the body where it is embedded) or the boundaries conditions (loads and displacements) only implies to modify the finite element model, without altering to any extent the current optimization-based design procedure. So, complicated, maybe three-dimensional, real industrial or engineering problems can be straightforwardly accounted for. For instance, to account for multiple loads implies an easy redefinition of the objective function. Let u 0 = u 0 (x, F (α) ) be the displacement caused by the external force F (α) in the domain Ω without the inclusion, and u = u(x, P, F (α) ) be the displacement caused by that force in Ω in presence of the inclusion Ω incl together with the cloaking device occupying Ω dev , where the material distribution is defined by P. In presence of multiple loads (2) ( ) load , the objective function representing the error in the accomplishment of the global cloaking task can be defined as the weighted sum of the RMSEs in the accomplishment of the individual tasks, that is where ω α is the weight assigned to the task accomplishment for the external force F (α) ; typically, ω α = 1/N load .
As an application example, let us consider the problem described in the previous section when a vertical tension is added, as shown in Fig. 4. Now, cloaking has to be achieved for horizontal compression (as in the previous example) as well as for vertical tension. In the holed plate fully made of nylon, i.e., without a cloaking device, the error in the accomplishment of the global task is = = .
= . The optimal device is that shown in Fig. 5 If we define a fully discrete device by assuming that those finite elements in the optimal device having an aluminum fraction greater than φ * are completely made of aluminum, the best design is that obtained for φ * = 0.033, shown in Fig. 5(j), for which the global error is f obj = 0.1101f nocloak , and the individual errors are = . RMSE 0 1120RMSE hcomp nocloak hcomp and = . RMSE 0 1081RMSE vtens nocloak vtens . If, for easier manufacturability, we eliminate from this device those little PTFE areas adjacent to the hole obtaining the device shown in Fig. 5(m), the accomplishment of the task is slightly affected, with f obj = 0.1104f nocloak , = . RMSE 0 1124RMSE hcomp nocloak hcomp , and = . RMSE 0 1084RMSE vtens nocloak vtens . Let us remark that, because of the highly non-discrete nature of the optimal device ( Fig. 5(g)), the accomplishment of the cloaking task is considerably affected by the black-and-white filtering. However, this is still satisfactory as it can be realized when computing the error Δ in the cloaking task as defined by Bückmann et al. 15 (see equation (10)): Δ = 10.74% for cloaking the horizontal compression, and Δ = 10.40% for cloaking the vertical tension.
Note that the displacements shown in Fig. 5 are those produced by the horizontal compression, say u hcomp . Since the devices (the optimal as well as the discrete ones) are practically symmetric with respect to the center, the displacements produced by the vertical tension, say u vtens , are easily derivable from u hcomp : actually, the horizontal and vertical components of u vtens have color maps almost identical to those shown in Fig. 5 but their thresholds are interchanged. Not surprisingly, the current devices accomplish the cloaking task almost equally for both loads.
It is also not surprising, and worth noting, that a device designed for multiple loads performs the cloaking task for each individual load poorly than the device specifically designed for such load, as done in the previous section.