Designing Shape Morphing Behavior through Local Programming of Mechanical Metamaterials

Shape morphing implicates that a specific condition leads to a morphing reaction. The material thus transforms from one shape to another in a predefined manner. In this paper, not only the target shape but rather the evolution of the material's shape as a function of the applied strain is programmed. To rationalize the design process, concepts from informatics (processing functions, for example, Poisson's ratio (PR) as function of strain: ν = f(ε) and if‐then‐else conditions) will be introduced. Three types of shape morphing behavior will be presented: (1) achieving a target shape by linearly increasing the amplitude of the shape, (2) filling up a target shape in linear steps, and (3) shifting a bulge through the material to a target position. In the first case, the shape is controlled by a geometric gradient within the material. The filling kind of behavior was implemented by logical operations. Moreover, programming moving hillocks (3) requires to implement a sinusoidal function εy = sin (εx) and an if‐then‐else statement into the unit cells combined with a global stiffness gradient. The three cases will be used to show how the combination of mechanical mechanisms as well as the related parameter distribution enable a programmable shape morphing behavior in an inverse design process.


Introduction
Often the shape of a device can be optimized for one specific functionality, but many applications serve conflicting functions. For example, airplane wings being efficient during cruising and is obtained by linearly increasing the amplitude of the lateral deformation. In the second case, the target shape is filled up in linear steps. In the last case, the target shape is achieved by shifting the position of a bulge in the material. To realize the first case, we need a unit cell that processes a linear function depending on one tunable geometrical parameter. For the second case, an if-then-else condition that can be tuned by a geometrical parameter was implemented in a unit cell. The third case requires a cell that processes a sinusoidal, nonmonotonic function in combination with an if-then-else condition as well as global stiffness gradient in the material. Figures 2 and 3 show the three unit cells that exhibit the required behavior (functions, if-then-else conditions) for the three aforementioned cases. They can be parametrized and assembled into metamaterials. Their complete geometrical description can be found in Section S1, Supporting Information. To achieve reversible behavior, we consider a purely elastic base material without plastic deformations. Global shape morphing requires to orchestrate the local deformation in a way that the target shape can be achieved when the metamaterial is deformed. One approach is to design the local PR that corresponds to a shape morphing behavior shown in Figure 1, case 1. For the sections of the part, where the PR is negative, the part contracts under compression and vice versa which can be designed in the unit cells. The hexagonal structure (Figure 2a) is well-known from literature [22] and has properties that strongly depend on the angle α. In the case of small deformations (ε << 1), the mechanical properties of such structures can be calculated analytically and are constant for a chosen angle. Such a representation was derived on the basis of classical homogenization concepts [23] (see Section S2.1, Supporting Information). The resulting linear stressstrain relation leads to constant stiffness and PR for small deformations. Figure 2b,c shows the derived Young's modulus E and PR ν for different angles α and a constant unit cell size. The variation of this angle allows to generate a broad range of properties only dependent on one geometric parameter. For α < 90°, the cell adopts the shape of a bow-tie and is auxetic. Choosing α > 90° leads to a honeycomb with positive PR. [15,16,18,19,24] Note that the adaption of the PR is accompanied by a change in the stiffness and vice versa. These properties can be decoupled through adapting other parameters such as beam thicknesses (t) or the base material.

Results and Discussion
Another possibility to design a target shape into the metamaterial is locally limiting the maximum deformation under a given load (see Figure 1, case 2). This can be achieved by abruptly increasing the local stiffness at a defined unit cell deformation. To gain these non-linear properties, we introduced a contact gap that can close and open during load cycles and that leads to a structural transformation (Figure 2d). The cell geometry changes between a bow-tie and one consisting of two trapezoids. This transformation can be considered in the homogenized model by introducing an internal variable δ that allows to distinguish between open and closed contact (see Section S2.1, Supporting Information). The linear stress-strain relation is now divided into two parts with different slopes. The constant properties of case a) establish an if-then-else condition:    b) and c) resulting constant stiffness and PR for different angles α over strain ε x . d) Geometry of the cell with inner contact. e) and f) Resulting stiffness for different gap widths μ g (μ g, 1 < μ g, 2 < μ g, 3 ) for fixed α. Under increasing strain, we see an abrupt change in the properties that can be interpreted as an if-then-else condition. which can be moved along the loading axis (see Figure 1, case 3), different mechanisms must be combined. Here, complex interactions between unit cells need to be programmed in the metamaterial. To design a suitable unit cell, the properties and geometric parameters in the former examples can be used as starting conditions. The cell must be extended so that large deformations are possible (see Figure 3a). In this case, the behavior cannot be derived analytically anymore but may be approximated with semi-analytical equations as well as numerically investigated using a finite element model (see Section S2.1, Supporting Information). We start from a honeycomb cell but change some geometric parts: We include hinges at the connection of the beams that allow them to rotate. The hinges are realized in terms of tapering of beams at the connection points. They are parameterized by the thickness t 2 that is smaller than the middle part of the beam t 2 < t. In consequence, the deformation is located in these parts of the beams. Additionally, a contact element has been implemented which is characterized by the angle β. It has a similar function as μ g in the upper example but blocks the rotation of the beams.
The angle α changes during uniaxial compression and the PR ν becomes a non-linear function (Equation 2) of the strain depending on the current angle (α 0 + Δφ) of the cell geometry (beam length l 2 , cell width 2b 0 , details can be found in Figure S1b, Supporting Information). Δφ describes the difference between α 0 (angle in the undeformed geometry) and α(ε x ) during loading. This leads to a non-monotonic lateral contraction visualized in Figure 3c. At low strains, the PR is positive ((α 0 + Δφ) > 90°) and changes to negative values for large strains ((α 0 + Δφ) < 90°). The deformation is then stopped through the contact that increases the stiffness and sets the PR approximately on zero. So three different ways of shape changing under load are combined in one cell: extension, zero deformation, and contraction. We can describe the cell's behavior as following: is based on a semi-analytical description assuming pure rotation of the beams which gives a rough estimation of the unit cell behavior (see Section S2.1, Supporting Information). Figure 3b,c shows the resulting stiffness and PR from the FEM simulation. The above-described parameters α, β, and t 2 can be used to control the processing functions. We can increase the stiffness by increasing t 2 , move the shift in the stiffness to lower (higher) strains by increasing (decreasing) β, and influence the start value of the PR through the parameter α. This allows to adapt single cells in an assembly. In a next step, macroscopic shape morphing structures were realized by creating 2D arrays of the unit cells shown in Figures 2 and 3 within a FEM framework. In the first structure, we varied the geometrical parameters α in a macroscopic material that consists of the hexagonal cells (Figure 2a). The corresponding data is displayed in Figure 4, first column. The second structure is composed of the cells shown in Figure 2d with a fixed angle α and a variation of the contact gap μ g (Figure 4, second column). For the third case, we show two examples that consist of the cell shown in Figure 3. The parameters α and β were fixed in both examples but t 2 varied linearly and sinusoidal over the x-coordinate of the samples (Figure 4, third and fourth column). Figure 4a shows the macroscopic scale with the applied boundary and load (details in Figure S3 Here, we plotted the y-strain over the x-coordinate which corresponds to the shape of the sample surface Γ* (marked with an orange line).
The parametrization allows us to control the behavior (linear, if-then-else, non-monotonic function), described in the unit cell section, locally and thus to program the global shape morphing. If we aim to achieve a specific shape morphing behavior in a material, we need to know how to adapt the geometric parameters in the individual cells by taking into account their respective stress/strain field. Due to the homogenized model of the hexagonal cell with and without contact, the ideal parameter distribution can be determined with optimization methods (for small deformations). To use this continuum descriptions for the materials, the size of the unit cell must be much smaller than that of the macroscopic specimen so that the macroscopic strain does not vary significantly over the length of a cell. However, comparison of the results for the fully resolved structure and the homogenized model reveals that, for the case in Figure 4a, good results can be achieved with as few as nine cells over the height of the structure, despite the fact that the angle varies strongly among the cells (see Figure S4, Supporting Information).
To gain a target shape, we aim to optimize, in particular minimize, the error or difference to a certain desired surface displacement (in the lateral direction) under compression. The error is also referred to as cost or objective function and can be specified as follows: where u and u are the actual and desired displacements, respectively. Here, u and u are functions of space, hence u ( ) J measures the error in the Euclidean norm. In order to minimize this error, we define γ as control, which is the design parameter as function of space. Therefore, we may write u = u(γ). For instance, we may consider the angle α, the contact gap μ g , or the hinge thickness t 2 in Figure 3 as control through the design. Note that in this work, we focus on only one control variable, more precisely we optimize w.r.t. to one design variable (e.g., either α or μ g , see Figures S4 and S5, Supporting Information). Nevertheless, optimizing with one design variable is still non-trivial, due to large number of unit cells in compound and, hence, large design space: Since γ is also a function of space, each cell (or finite element) is equipped with a degree of freedom within the parameter space. Therefore, optimizing γ within a compound domain means the simultaneous optimization of the design parameter in each unit cell. An adjoint sensitivity approach is a convenient choice for solving such optimization task. Since we want to apply an iterative gradient descent method, [25] the functional or Fréchet derivative ( ( )) u J γ δγ ′ is necessary. Roughly spoken, it may be understood as directional derivative with δγ as variational or virtual term. In order to shorten the notation, we use the abbreviation In this work, we apply the adjoint sensitivity approach; that is, before each gradient descent iteration, the original mechanical balance of equilibrium and its adjoint version have to be solved. For these two steps, we use the finite element software CalculiX. [26] The resulting solutions are then processed for the evaluation of the total parameter derivative d γ J . Roughly, one can describe an iteration like this: The iteration above is repeated until convergence occurs (Details in Section S3, Supporting Information).
In our first example, the honeycomb cell was used and the parameter α was varied in the range of 60° to 120° to have positive PRs in some regions as well as negative in others. Hence, the lateral contraction is non-uniform over the x-coordinate of the sample but its amplitude is linear in the applied horizontal global displacement u x (Figure 4d, left). As we optimized the parameters, we could reduce the range to values between α = 84° -94° (Figure 4c, left). Ref. [18] as well as [19] report 3D-printing and analysis of similar structures (but for large deformations) resulting in similar shapes as shown in Figure 4d, left. However, we choose a multiscale approach which allows us to treat the structures as materials (with a lot of unit cells) and to transfer the methodology to more complex problems (e.g., include more parameters). This case shows how a predefined target shape, with an amplitude depending on the load, can be achieved through an adaption of geometrical parameters.
The second structure (Figure 4b  gap μ g was locally adapted in the following. The achieved lateral strain (ε y ) over the x-coordinate of the part looks similar to the first example, but results in a non-uniform (in this case negative) lateral contraction. If we try to relate the lateral contraction to the load amplitude, we do not get a linear relation anymore. Figure 4d, left, shows different shapes for different applied global displacements u x . This is a result of the local ifthen-else conditions implemented in the unit cells. The local stiffness change leads to roughly one order of magnitude stiffer regions. They deform less under increasing amplitude of the global displacement and in addition have a positive PR. From then on, these cells are (locally) stress-controlled and their shape remains practically unchanged (but the others cells can still be deformed). In consequence, the surface displacement is a piecewise-linear (non-smooth but continuous) function of the global displacement. Every time a unit cell changes its contact state, the global shape changes as the smooth deformation is blocked. When all contacts are closed, the overall stiffness of the material increases significantly.
For this metamaterial, it is possible to optimize both the surface displacement and the overall stiffness but also to control a deformation path. In Figure 4d, middle left, we show how the shape changes from flat to non-smooth, to a "u-shape." These structures are interesting for applications in which stiffness as well as shape should be adaptable (friction change, packaging, etc.). The optimization of μ g leads to a narrow parameter distribution in the material (Figure 4c, left), wherein the unit cells must be adapted only slightly. Intended stiffness changes can be found in few metamaterial structures. Fang et al. show an overall stiffness change in a self-locking origami structure [27] and in ref. [28] externally controlled electromagnets have been used to open and close structures. Another possibility is to exploit instabilities to have a structural transition. [29,30] Moreover, bistable structures can also lead to such changes combined with memory behavior as shown in ref. [31] and ref. [32]. These approaches lead to noncontinuous functions due to the snap-through. Also in granular system, we can introduce stiffness changes due to particle jamming. [33,34] However, our structure reacts passively under strain control, is simple to realize, and allows continuous deformations as well as a local implementation. We only design a jump in the local properties (ν, E) but not in the displacements. So we avoid instable systems which are complicated to handle, but could be useful to change abruptly between two different shapes. In this example, the shape is not only programmed by the choice of the initial geometry but also by the loading amplitude.
In the third case, the aim is to program a mechanical metamaterial to form a surface hillock similar to the first example, but now, this hillock should move along the x-axis depending on the applied strain in that direction. This requires to implement a non-monotonic behavior and to consider large deformations. For this reason, we used the unit cell shown in Figure 3 that can expand as well as contract and assembled an array of 2×10 cells. The hinge thickness was defined in each cell so that the stiffness increases from left to right (t 2 = t 2,0 + t 2,0 x, Figure 4c, middle right). The cells are connected in series (with increasing local stiffnesses) and in parallel (with equal local stiffnesses). The non-monotonic, strain-dependent stiffness of the unit cell leads to a global stiffness minimum which moves through the material. The deformation of the structure starts on the left and migrates to the right when the stiffness of the left part increases. The non-uniform local strain leads to different lateral contractions in the cells depending on t 2 . Due to the local non-monotonic behavior of the single cell, a nonlinear, non-monotonic behavior of the system can be achieved. Figure 4d, middle right, shows the different surface morphologies, which can be controlled with this structure. If we choose the hinge thickness according to a sinusoidal distribution (t 2 = t 2,0 + t 2,0 sin (xπ), Figure 4c, right), the surface can even have two bulges that merge together after a certain strain (Figure 4d, right). The variation ranges from wedged and warped to a plane surface shape although the design space is not widely explored yet. Videos of the FEM simulations of the two structures are provided as Videos S7 and S8, Supporting Information. In the future, our homogenization and optimization framework shall be expanded to include large deformations, history dependency, and more geometric parameters to handle this complex behavior. Just like in the second example, the initial geometry as well as the amplitude of the global displacement can be used in this structures to program the shape. Moreover, the more complex unit cell geometry allows a larger range of target shapes possible within one sample. Producing such kind of materials is challenging due to their complex, inhomogeneous geometries and the difference in size between device volume and unit cell ( unit cell volume device volume << 1). A first approach was made by 3D-printing structures of a few unit cells, but there are still a lot of limitations like nozzle sizes and overall build volume. Nevertheless, we produced a demonstrator part made of the thermoplastic elastomer (TPE) Elastollan C78 A15 corresponding to the structure shown in Figure 4, third column, with a linear distribution of the parameter t 2 (see Section S4, Supporting Information). As an experiment, the demonstrator has been fixed on a plate so that a rolling contact in x-and y-direction occurs. This permits a qualitative validation of the simulation results. In Figure 5a,b, finite element simulations (with ABAQUS 2018, [35] Section S2.2, Supporting Information) are shown next to the 3D-printed part under different load conditions. Figure 5c shows the comparison of the shapes of the printed samples and FEM simulations. Here, we used the local lateral strain over the x-coordinate of the sample under different global displacements (u x,1 < u x,2 ...) to describe the deformation and shape of the sample (see Video S7, Supporting Information). The experimental data was obtained by optical evaluation. To gain a better quantitative accordance complexer material models, the tolerances in the 3D-printing process as well as the non-ideal boundary conditions in the experimental evaluation have to be considered.

Conclusion
All in all, mechanical metamaterials open up a multiscale design space which can only be utilized if the functionality can be related to their inner structure and vice versa. We have shown that programming elements can be used in the design process of materials to obtain predefined behavior. Designing programmable behavior into such metamaterials requires to introduce logical operations (e.g., if-then-else), as well as the capability of processing functions (e.g., ε y = f(ε x )), which can be parametrized. The selection of such elements as well as the optimization of parameter distributions is required to integrate certain (shape morphing) functions into materials. Besides controlling the final shape of a deformed material, the course and intermediate shapes can be programmed as shown in cases 2 and 3. With increasing number of controllable parameters within the unit cells, complex, non-linear system behavior can be implemented. This level of control allowed us to create a material that behaves similar to the locomotion of earthworms (see Figure 5d) which is known for its robust and space-saving way to move. [36] A locally controlled shape morphing and a material that can expand as well as contract is required. The different parts in the material have to expand and contract alternately. In the field of soft-robotics, this peristaltic is controlled externally through, for example, magnets [37] or motors [5,36] and results in complex technical solutions. With our approach and the material shown in Figure 5, we can transfer this functionality directly in the material and control it through the input load. We also see potential in using such materials to control transport processes on surfaces or in tubes where purely mechanical solutions would also lead to a simplification of the systems. In our example, we varied only one of the several parameters in the unit cell; however, we see potential of the exploration of multiple parameters. This could -open doors to control shapes in aerodynamic applications within the material as well as designing adapted fixing systems such as plugs and pins that could also profit from local shape control. Nevertheless, we have to integrate large deformations in our framework and to deal with straindependent target functions. Therefore, we need to develop a data-based approach wherein all representative load cases will be pre-calculated. Furthermore, the design process needs to be optimized to implement logical elements in unit cells while considering manufacturing restrictions and avoiding instabilities. Nevertheless, such instabilities can also be used as design elements. [29,38] Overall, our multiscale approach to use parametrized unit cells, homogenization techniques, and optimization leads to a reduction of the complexity and helps to identify programmable materials to replace complicated technical systems.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.