Computational design of 3D-printed active lattice structures for reversible shape morphing

Abstract

Active structures can adapt to varying environmental conditions and functional requirements by changing their shapes and properties, which makes them suitable for applications in changing environments as found in aerospace and automotive. Of special interest are light and stiff structures with shape morphing capabilities, which is naturally contradictory. Existing concepts in literature can be limited to a single, non-reversible actuation and are difficult to design due to the inherent complexity of large-scale lattices with many elements and complex target deformations. Here, we show how 3D-printed active materials can be combined with an efficient computational framework to design large-scale lattice structures that can change their shape between an initial state and a target state. The reversible deformation is controlled by a single actuation input and heating of the structure. Numerical and experimental results show the generality of the proposed method and the applicability to different problems such as morphing airfoils.

Graphic Abstract

Introduction

The concept of shape morphing is abundantly found in nature, from adaptive leaves of plants [1] to morphing wings of insects and birds [2, 3], and likewise imitated by engineers and scientists to create structures that can adapt to different operating conditions. In general, these active structures react to changing environmental conditions and change their shape and properties to meet new functional requirements. This is often enabled using active materials, which can change their properties in response to external stimuli such as heat, moisture, and light [4]. The recent popularity of active materials in engineering and science is partially owed to the possibility of 3D printing these materials, combining their exceptional properties with the inherent geometric and topological complexity of advanced 3D-printed structures. Among others, this combination constitutes the field of 4D printing [5], where temporally changing objects and structures are explored that can fulfil different tasks such as autonomous cargo delivery [6], self-assembly [7, 8], and locomotion [9]. In contrast to conventional morphing structures, which often require complex assemblies of many individual parts, 4D-printed structures can be fabricated monolithically where functional components such as compliant hinges and active elements are directly integrated. This reduces the overall number of parts, assembly time and maintenance requirements, as, for example, no lubrication is required for compliant hinges. This makes active structures especially interesting to industries where these factors are important, such as aerospace and automotive [10, 11]. However, the development of these structures is still limited by a lack of efficient computational design methods, since many examples found in literature are designed by hand and thus are limited by engineers’ imagination and experience. Additionally, shape morphing structures are often restricted to planar 2D structures or are made of compliant materials that do not possess any advanced load-carrying capabilities. To make them attractive for real-life applications, for example in morphing airfoils and car panels, active, shape-morphing structures must be simultaneously rigid, load-carrying and compliant, which is naturally contradictory. Finally, some 4D-printed structures are designed for one-way usage and “re-programming” after actuation can be difficult, especially for structures with many individual actuators and complex deformation patterns.

Here, we combine the advantages of 3D-printed active materials and advanced computational design and modelling to design lightweight lattice structures that can morph between an initial state and a target state in a reversible way. This is achieved by exploiting the thermomechanical properties of two different materials, which have the same stiffness at room temperature but different stiffnesses at higher temperatures. By optimizing the material distribution in a predefined ground structure, stiffness gradients emerge at higher temperatures and enable global shape morphing, which can be controlled by a single input actuation displacement. All deformations are based on the reversible shape-memory effect of the materials, which leaves the initial topology always intact and ensures that the structures remain mechanically sound.

Background

Materials

The basis of many active structures that are activated through their materials is a combination of multiple active materials with different mechanical and thermal properties, where for example the stiffness and strength vary over several orders of magnitude [4]. Modern additive manufacturing processes enable the simultaneous deposition of multiple materials at the microscale, which enables the fabrication of composite materials with tailored properties [9, 12]. Similar to classic composite materials, the interaction between the different components can yield unique properties that cannot be achieved by the constituent materials individually. A class of materials that is frequently found in literature on active structures are shape-memory polymers (SMPs) [13]. SMPs are available for many 3D printing processes such as material jetting [14], stereolithography [15], and fused deposition modeling [16]. The mechanical properties of SMPs are directly linked to their respective glass transition temperature (T_g). Below the T_g, the polymer is in its “glassy” state, where it is normally relatively stiff and brittle. Above the T_g, the material is in a viscous or “rubbery” state and the stiffness of the material decreases. In practice, the transition temperature is mostly not a discrete temperature but a temperature range in which the properties of the materials change. When a SMP is deformed in its rubbery state, the deformations can be fixed by cooling the material below the T_g. This process is often called “programming” as the material is programmed from one initial shape into a new shape. In contrast to many other materials, SMPs have the ability to remember the initial shape and recover it on re-heating above the T_g. The deformations and strains during this reversible programming and recovery cycle can be relatively large, but in general do not damage the material itself. Hence, the shape memory effect can be used to store and release deformation energy in the materials, which then can be used as actuators or active elements in shape-morphing structures [14]. While traditional SMPs require an initial programming step before they can be actuated, strategies exist to overcome the limitations of the one-way actuation. This involves polymer laminates of multiple materials [17, 18], exploiting multi-shape memory effects of materials with a large transition temperature range [19], and combining materials with different stimuli such as heat and moisture [20, 21]. While all these strategies involve additional materials and geometric control, a true reversible shape memory effect is found for example in liquid crystal elastomers (LCEs) [22]. However, as these structures are currently still very expensive, difficult to fabricate, and the shape memory effect can be unstable [23], they are unsuitable for the fabrication of active structures at a larger scale, as shown in this work.

Shape-morphing structures

The usage of different materials with distinct geometrical designs enables the actuation at different temperatures for sequential actuation, which can be used to design multi-state structures [24] and autonomously manoeuvring systems [6]. The drawback of these structures is, however, that the actuators have to be programmed manually before actuation, which makes the process non-reversible and infeasible for more complex structures with hundreds or thousands of active parts. Such structures require the support of computational tools, which can help to design structures with complex deformation patterns and shapes at a large scale. However, active structures rely on the interrelation of micromechanical and macromechanical effects, which have to be modelled with sufficient accuracy to generate viable results. This can involve complex nonlinear material models [25] or prestresses [26], which makes it computationally expensive to solve larger problems. More studies focus solely on the optimization of shape-morphing structures without active materials. Common techniques involve the optimization of geometric parameters of lattices [27], the mapping of locally deformed, precomputed patterns to achieve global deformations [28] or multi-stage approaches where topology optimization and mapping strategies are combined [29]. However, these shape transformations are purely elastic and cannot be fixed permanently. Weeger et al. propose an approach where lightweight lattice structures and active materials are combined to design structures that can be programmed into complex 3D shapes [30]. Even though the target shapes can be fixed and the programmed deformations are reversible, the structures are, in general, compliant in both states due to the use of both rigid and compliant materials, which makes this combination less feasible for applications where load-carrying structures are required.

In this work, we combine the advantages of selected SMPs and their distinct thermo-mechanical properties, lattice structures, and advanced computational methods to efficiently design 3D-printed structures that can reversibly change their shape from an initial configuration to a second target shape by a simple input actuation while remaining structurally intact. We start by investigating the thermo-mechanical properties of the materials used in this work, from which we derive the actuation concept of the structures. We then derive an optimization problem formulation and show how it can be applied to different example structures that are 3D printed and tested experimentally.

Results

Shape morphing cycle

In contrast to shape morphing structures where individual elements or bars have to be actuated separately to achieve a predefined target shape, we introduce a concept where uniform heating activates the structure and enables shape morphing. The concept is based on exploiting the different mechanical properties at different temperatures of materials that are 3D printed with the Objet500 Connex3 3D printer (Stratasys, Eden Prairie, MN, USA). The 3D printer uses the material jetting technology and allows for printing multiple materials at the same time. Figure 1a shows the Young’s modulus at different temperatures of the two materials “VeroWhite+” (VW) and “High Temperature” (HT). The values are obtained by tensile tests according to the respective ASTM norm [31], conducted on an Instron ElectroPuls E3000 mechanical testing machine and an Instron 3119 environmental chamber.

Figure 1

(a) Temperature dependency of the Young’s modulus of the two materials “VeroWhite+” and “High Temperature”. (b) A small 2D lattice structure with compliant members (gray) and rigid members (black). Depending on how the rigid members are arranged in the structure, the structure can have negative Poisson’s ratio (solid black lines) or positive Poisson’s ratio (dashed black lines) when expanded horizontally. (c) The structure from (b) in the deformed state under horizontal expansion for both material distributions. (d) Schematic of a larger lattice structure (gray) with supports (left) and input actuation displacement (right). A tailored material distribution inside of the lattice structure can enable complex target deformations such as the sine shape at the top and at the bottom. (adapted from Ref. [32]).

The markers indicate the mean value of five samples, and vertical error bars indicate the standard deviation. At room temperature and up to about 50 °C, both materials VW and HT have approximately the same Young’s modulus and decrease almost linearly. At higher temperatures, the stiffness of VW drops, while the stiffness of HT further decreases, almost linearly. Even though both materials become relatively compliant, the stiffness of HT remains higher than the stiffness of VW at temperatures up to 80 °C. This behavior is caused by the two different glass transition temperatures of about 55–60 °C (VW) and about 70 °C (HT).

Instead of achieving a target deformation by recovering a pre-programmed state, as often seen in literature, we use this difference in the thermo-mechanical properties of the materials to introduce local stiffness gradients in the structures at higher temperatures, which cause the target deformations under a single actuation input displacement. A similar concept was introduced by Weeger et al. for semi-compliant, voxelized beam structures [30]. Figure 1b schematically shows a 2D lattice structure where each bar either consists of VW (gray lines) or HT (black lines). The solid black line and the dashed black line show two different possible distributions of HT material. At room temperature, all bars have the same Young’s modulus of about 2100 MPa and the lattice has a homogeneous stiffness distribution. At 80 °C, the stiffness of the HT bars (black) is reduced to about 305 MPa and the stiffness of the VW bars (gray) is reduced to about 8 MPa. Since the HT bars are about 38 times stiffer than the VW bars, local stiffness gradients are introduced to the structure and the overall mechanical response is now mostly governed by the distribution of the stiffer HT material. Hence, by controlling the material distribution in the structure, the mechanical properties and the deformation under external displacements can be controlled.

In the first case, HT material is deposited as indicated by the solid black lines and forms a re-entrant cell, which is generally associated with auxetic behavior [33]. Hence, at high temperatures a horizontal expansion of the structure will result in a vertical expansion (Fig. 1c top), the overall Poisson’s ratio of the structure is negative $\nu < 0$. In the second case, HT material is deposited as indicated by the dashed black lines and forms a diamond cell, which will contract vertically under a horizontal expansion with $\nu > 0$ (Fig. 1c bottom). This simple example shows that it is possible to control the overall deformation of a structure at high temperatures by assigning the two materials VW and HT to the individual bars. Figure 1d shows schematically how this concept can be implemented in a mechanical structure. The gray region represents a general lattice structure, which is fixed at the left side and which can be mechanically deformed by a local input displacement at the right side, for example via a pneumatic actuator. The dash-dotted lines indicate the global deformations of the structure at high temperatures, which are determined by the exact material distribution in the base lattice.

By controlling the input deformation and the temperature of the structure, this setup enables reversible shape morphing behavior. Figure 2 schematically shows a 2D lattice structure and the respective stiffness distribution during all states of the shape morphing cycle. In step one, the structure is at room temperature (23 °C) and all bars have the same stiffness, independent of the material distribution. In step two, the structure is heated up to 80 °C. The bars that consist of VW, indicated by the gray lines, become more compliant than the bars that consist of HT, indicated by the black lines.

Figure 2

All steps of the shape morphing cycle for the “2D Sine” structure (adapted from Ref. [32]). All bars have the same Young’s modulus at room temperature (1). At 80 °C, the stiffness of the members changes due to the two different materials used (2). An input displacement is applied to the right side of the structure and the distinct material distribution causes a predetermined global deformation, for example a sine shape here (3). The structure is cooled down, which fixes the deformations. (4) The initial shape can be recovered by re-heating the structure due to the shape memory properties of the materials.

Next, an input displacement is applied locally to the right side of the structure. Due to the heterogeneous stiffness distribution, the structure deforms into the desired shape, which is the sine shape at the top and at the bottom here. To fix the deformations, the structure is cooled down to room temperature in step 4 and the stiffness in all bars is again the same. Since the overall deformations are based on deformations in the members, the overall topology of the structure is preserved. Due to the shape memory properties of the materials, the applied and fixed strains in the bars are reversible, and the original shape can be recovered when the structure is heated once again to 80 °C.

Since different material distributions yield different deformations for given input displacements, it can be difficult to find a material distribution for which the structure deforms to meet the target shape. For small lattices, intuitive configurations such as the re-entrant or diamond configurations in Fig. 1b might yield qualitatively good solutions. However, a computational design framework is required to systematically design larger scale lattices with more complex deformations while taking into account that the optimization problem is nonlinear and non-convex. In this work, we develop a two-step computational design optimization framework to solve the discrete material distribution problem. To make use of gradient-based optimization techniques for efficiently solving large-scale problems, the discrete problem is relaxed and solved with an interior point line search filter method (IPOPT) as a first step using a random starting point. As the relaxed problem can yield solutions that are not fully discrete, i.e. some bars with intermediate material properties remain, in a second postprocessing step, a Genetic Algorithm (GA) is used to find fully discrete solutions based on the almost discrete solutions from the first step. The optimization model and methods are described in detail in “Methodology”.

Optimized Structures

To show the applicability to complex 3D problems and the generality of the method, Fig. 3 introduces the three example structures “2D Sine”, “3D Sine”, and “Airfoil”. The example “2D Sine” in Fig. 3a consists of 6 × 4 square unit cells with 154 bars, where each unit cell has the side length $L = 17.0\;{\text{mm}}$ and a central node in the middle. Hence, the length of the horizontal and vertical bars is $l = L$ diagonal bars is $l = (\sqrt 2 /2) \times L$. The $n_{{\text{T}}}$ target displacements at the top and at the bottom, indicated by the dashed blue lines, are defined on in total 10 individual nodes and describe mirrored sine curves with an amplitude of a = 3.0 mm. The input actuation is a displacement of all nodes at the right side of the structure, indicated by the black arrows. The structure “3D Sine” consists of 4 × 4 × 2 octet-truss unit cells with 896 individual bars. All bars are of equal length $l = 14.1\;{\text{mm}}$. In contrast to the 2D example, this structure is deformed by compressive displacements at the top. The target displacements on the left side and on the right side describe mirrored sine shapes with an amplitude of a = 3.0 mm and are defined on in total 46 individual nodes. The structure “Airfoil” has the shape of the back half of a NACA 6420 airfoil [34]. This shape is extruded in z-direction and discretized with 18 × 2 × 8 body-centered cubic unit cells, which results in 3516 bars in total. Due to the mapping of the unit cells to the airfoil geometry, the bar lengths $l$ vary between 12.4 mm at the root and 0.7 mm at the tip. The target shape describes a deflection of the tip of the airfoil downwards by a = 10.0 mm, as it could be used for example for shape-morphing wings. The target shape is defined by, in total, 105 individual nodes at the top and at the bottom of the airfoil. A bar radius of r = 0.75 mm is chosen for all structures.

Figure 3

Lattices and optimization problem setup of three example structures with boundary conditions, input displacements (arrows), and target displacements (blue dashed lines/blue surfaces). (a) “2D Sine”, (b) “3D Sine”, (c) “Airfoil”.

Each example structure is optimized 10 times, due to the random starting point, and the arithmetic mean and the standard deviation for both optimization steps are given in Table 1. For all three example structures, the best solution after both steps is selected out of the 10 optimization runs, fabricated and experimentally tested. It can be observed that the solution of the relaxed problem after the first optimization step achieves lower objective function values for all three examples than after the postprocessing step, as the restriction on purely discrete values in the postprocessing step further constrains the design space. On the left side of Fig. 4, the best optimized solution for each example is shown in the deformed state. The deformations are computed by a linear-elastic, truss FE analysis with material properties at 80 °C. The colors of the bars indicate the material distribution, where light gray represents VW material and dark gray indicates HT material. The transparent blue dashed lines and the transparent blue planes indicate the predefined target displacements, respectively. In general, the target shapes are matched by the optimized structures well. However, the absolute values of the objective function of the three examples after the postprocessing step are significantly different, as summarized in Table 1.

TABLE 1 Statistical results of ten consecutive optimization runs for each example.

Figure 4

Optimized structures in the deformed state, as predicted by the FE analysis (left) and fabricated structures (right). Black members represent the stiffer HT material and gray members represent the more compliant VW material. Blue dashed lines and blue surfaces indicate the target shapes. The structures are the solutions with the best objective function value f after 10 consecutive optimization runs, respectively. Further, the optimization mean error $e_{{{\text{opt}},{\text{mean}}}}$ and experimental mean error $e_{{{\text{exp}},{\text{mean}}}}$ are shown. (a) “2D Sine”. Blue lines in the right image show the numerically predicted deformed shape of all bars. (b) “3D Sine”, (c) “Airfoil”. The solid green line indicates the undeformed shape (scale bar 1 cm).

This stems from the definition of the objective function $f$, i.e. the sum of the squared error over all target nodes. Larger structures with more target nodes such as the “Airfoil” more likely have larger objective function values than smaller structures with fewer target nodes. To make these values comparable, we compute the “mean error” of the optimization at every node via $e_{{\text{opt,mean}}} = \sqrt {f/n_{{\text{T}}} }$, with the best objective function value after the postprocessing step $f$ and the number of target displacements $n_{{\text{T}}}$. This metric can be interpreted as the mean error at every target node, i.e. the mean distance between the deformed shape predicted by the simulation and the target shape, evaluated at all target nodes of the optimized structure. The corresponding values are found in Table 1 and show that the optimization mean error at all nodes is between $0.12$ and $0.39\,{\text{mm}}$, which is the same order of magnitude for all examples.

To experimentally validate the feasibility of our approach, the selected structures are 3D printed on a Stratasys Objet500 Connex3 printer with the default options of the digital material mode. The structures “2D Sine” and “3D Sine” are heated up to 80 °C in an Instron 3119 environmental chamber. The actuation displacement is applied to the structure “2D Sine” by an Instron ElectroPuls E3000 mechanical testing machine. Solid material is added to the left and to the right of the structure, where the tensile grippers of the testing machine are attached. The structure “3D Sine” is deformed in the same testing machine using compression plates. Since the structure “Airfoil” is too big to fit in the environmental chamber, the structure is mounted to a rig, immersed in a 400 mm × 400 mm × 400 mm glass container filled with 80 °C hot water, and actuated manually with a metal rod. On the right side of Fig. 4, the fabricated structures are shown in the deformed state at 80 °C. Both fabricated materials VW and HT have the same color and cannot be distinguished in the 3D-printed structures. The target shapes are indicated by transparent blue dashed lines, respectively. In Fig. 4a, the full deformed structure as computed by the FE analysis is additionally shown on top of the fabricated structure with transparent blue lines. For better comparability, the fabricated structures in Fig. 4b and c are shown in the front view. In Fig. 4c, the thin solid green line indicates the initial shape of the “Airfoil” structure.

The deviations between the experimentally deformed structures and the numerically predicted shapes at the target nodes are determined with the image processing software ImageJ [35]. For the 3D structures, the experimental error is measured at all visible nodes in the front view. The mean error at the target nodes $e_{{{\text{exp}},{\text{mean}}}}$ is reported in Fig. 4 on the right side below the images and in Tab. 1. Even though the simulated structures match the target shapes with almost no visible deviations, small deviations on the millimeter scale are visible between the target shapes and the fabricated structures. The mean error is smaller than 0.6 mm for the examples “3D Sine” and “Airfoil”, but it is significantly larger for the structure “2D Sine” with 1.28 mm. Further it can be observed that individual beams buckle, which does not occur in the linear-elastic truss FE simulation.

The material distribution in the structure “2D Sine” is almost symmetric with respect to the horizontal symmetry line. The more rigid HT material, which dominates the structural behavior, forms a re-entrant shaped cell in the left part of the structure and causes the structure to expand vertically due to its auxetic behavior. In the right part of the structure, HT material is found at the top and bottom, which are connected to the actuated nodes on the right side by diagonal lines of HT material. Like this, these regions of the structure are “pulled” towards the center and the structure contracts. Between the left side and the right side of the structure, a vertical line of HT material prevents large vertical displacements of the top and bottom. The material distributions in Fig. 4b and c are more complex and geometric patterns that are related to the deformation behavior of the structures as in Fig. 4a are not directly visible anymore.

Figure 5 shows the convergence behavior of the three optimization runs during both optimization steps of the structures in Fig. 4 and the corresponding computation times. All computations are carried out on a commercially available Dell XPS 15 9500 notebook with an Intel Core i7-10750H processor. The computation times are in the order of magnitude of seconds to minutes, and the computation time increases with the complexity of the problem. Oscillations at low objective function values, as seen for example in the black “Airfoil” line, occur as the algorithm tries to find mechanically valid configurations while simultaneously reducing the objective function. The objective function values increase between the end of the first and the beginning of the second step as the remaining intermediate values are discretized. This effect is most noticeable in the “Airfoil” example as its objective function is generally higher than for the other examples.

Figure 5

Convergence behavior of the three examples “2D Sine”, “3D Sine”, and “Airfoil” shown in Fig. 4 during both optimization steps, and the corresponding computation times.

To experimentally show the reversibility of the programmed deformations, the structure “3D Sine” is cooled down for about 1 minute after its deformation at 80 °C by removing the environmental chamber. Next, the upper compression plate, which was used to apply the input displacements, is moved back up to the initial position. The structure is now in a deformed, stable state where both materials have again the same Young’s modulus of about 2100 MPa. To recover the initial shape, the environmental chamber is reinstalled, heated up to 80 °C, and the recovery of the structure is visually recorded. Figure 6 shows the recovery of the “3D sine” structure at 80 °C. The structure is initially fixed in the deformed shape and stable at room temperature. Upon heating up to 80 °C, the structure slowly recovers its initial shape, which is reached after about 6 min.

Figure 6

Recovery of the initial shape of the “3D Sine” structure by re-heating it up to 80 °C (scale bar 1 cm).

To investigate the mechanical behavior of structures in the deformed state, the structure “3D Sine” is mechanically tested at room temperature in both the initial state and in the deformed state, respectively. One structure in each state is compressed vertically until failure in the Instron ElectroPuls E3000 testing machine with a testing speed of 2.5 mm/min. The load and displacement data are recorded. Figure 7 shows the load carrying capacities of the “3D Sine” structure in the initial state (solid black line) and in the deformed state (dashed black line) at room temperature. In its initial state, the structure can carry a maximum compressive load of about 1250 N, while in the deformed state the maximum load is reduced by about 62% to about 470 N. The maximum displacement before failure however increases from about 2.6 mm to about 4.4 mm. While the load-displacement curve is theoretically linear at the beginning, some slack is observed here as not all parts at the top and bottom of the lattice are immediately touching the compression plates of the testing machine due to small fabrication inaccuracies. The blue solid and blue dashed lines indicate the stiffness of the structure in the initial and in the deformed state as predicted by the linear-elastic, truss FE simulation. It can be observed that the stiffness is almost the same, which indicates that the reduction in stiffness observed in the experiments does not stem from the shape change but rather from the buckling of the individual struts.

Figure 7

Experimentally determined mechanical behavior of the “3D Sine” structure in the initial state (solid black line) and in the deformed state (dashed black line), as well as the linear-elastic stiffness in the initial state (solid blue line) and in the deformed state (dashed blue line) as predicted by the FE simulation.

Discussion

The results in Fig. 4 and Table 1 show that the proposed concept for reversible shape morphing and the optimization framework can accurately produce structures that deform according to the defined target displacements and boundary conditions. With mean errors smaller than 0.4 mm, the deviation between the simulated, optimized structure and the target shape is orders of magnitude smaller than the size of the structures. Further, the mean errors between the fabricated structures and the target shape, which are between 0.36 and 1.28 mm, are also orders of magnitude lower than the overall structural size. In general, this increase of the error between optimization and experiment shows that the fabrication and testing inaccuracies have a larger influence on the final shape of the fabricated structures than the computational error. This includes mostly printing inaccuracies such as micro-defects, imperfections in applying the boundary conditions such as non-homogeneous clamping, and inaccuracies in setting and measuring the environmental temperature. The especially high error values of the "2D Sine", however, can mostly be explained by the specifics of the 2D model. As the structure in the experiment is not supported in the z-direction, it can warp out of plane. This out-of-plane warping is observed during the experiment and reduces the in-plane deformations, respectively. As such, the largest errors are seen in the "2D Sine" example in the middle and in the left part of the structure, as the struts that are supposed to move "towards the outside" warp out of plane. Without any support in the z-direction, out-of-plane warping can be caused already by small fabrication inaccuracies or inaccuracies in setting up the experiment. As this is not a problem in 3D structures, the errors are significantly lower. Hence, improving the fabrication and testing methods would have the biggest effect on improving the overall accuracy of the approach right now, which is, however, limited by the currently available technologies to print multiple active materials simultaneously.

The example “2D Sine” shows that the algorithm can find solutions with patterns that resemble shapes such as the re-entrant cell, which are intuitively expected in regions with negative Poisson’s ratio. The symmetry of the problem is reflected by the symmetry of the solution, even though a fully symmetric solution is not obtained. With increasing complexity of the examples “3D Sine” and “Airfoil”, no distinct patterns and no symmetric material distributions are directly visible anymore. This behavior is expected due to randomness introduced by the random starting points and the stochastic postprocessing step. Hence, finding a symmetric global optimum is not guaranteed.

The influence of the starting point on the quality of the solution can be estimated by comparing the mean and the SD of the objective function during the IPOPT optimization in Table 1. For the examples “2D Sine” and “3D Sine”, the optimization yields repeatedly “good” solutions, indicated by the relatively low SD. The influence of the starting point is small. The relatively large SD of the “Airfoil” example indicates that the quality of the solutions strongly varies for the ten optimization runs. Hence, the starting point strongly affects the results and must be taken into account to generate “good” designs, for example by running multiple optimizations.

The two-step optimization approach allows for efficiently solving problems with large numbers of variables. Figure 5 shows that IPOPT converges to a valid solution for all three examples in under 3 min computing time, and the postprocessing step takes less than 5 min for the most complex “Airfoil” example. Since the postprocessing step is computationally expensive compared to the IPOPT optimization step, we evaluate its usefulness by solving the “Airfoil” example as the largest example shown in this work with the stochastic GA only, i.e. from a random starting point without the first IPOPT step. Table 2 shows the mean objective function of 10 optimization runs and the mean computation times of the postprocessing step and the solution obtained by only using the GA. The mean objective function of using only the GA is more than three times worse than the mean objective function of the postprocessing step. Since the same settings are used, the computation time per iteration is about the same. The overall mean computation time of only using the GA however is about 1.5 min lower. This shows that the GA converges quicker to a local minimum, which, however, is worse than the solutions found with a better starting point. To show the computational cost of finding solutions similar to the postprocessing step by only using the GA, we increase the population size from 150 to 500. With this setting, the mean objective function value can be reduced to 55.91, which is similar to the postprocessing step, but the mean computational time increases to more than 16 min. Hence, while using only stochastic optimization methods, as also shown in literature [32], can be feasible for smaller problems such as the 2D sine, gradient-based methods or a combination of both as proposed in this work increase the computational efficiency and the accuracy of the solutions for larger problems with many variables.

TABLE 2 Mean objective function value and mean computation time of ten optimization runs for the “Airfoil” example of the postprocessing step.

The FE discretization, which the optimization is based on, models the bars of the structures as pin-jointed trusses. This simplification assumes that loads are only transmitted as tensile and compressive loads and no bending moments occur. As the 3D-printed structures are printed without pin joints, the bars are fused together at the nodes and the bars can deform in bending. This behavior would be in general better captured by beam models. However, by limiting the choice of ground structures to stretch-dominated ground structures, the truss model can still be accurate. Stretch-dominated structures are mostly loaded in tension and compression, i.e. the individual bars do not have to transfer any bending moments due to the redundancy of bars. Hence, a truss model can describe the structures with sufficient accuracy, which is confirmed in this case by the overall good agreement between the simulation and the experiments. Another characteristic of stretch-dominated structures is that they are statically indeterminate, i.e. they have redundant bars [36]. Choosing an overdetermined topology gives the optimization algorithm the possibility to form “substructures” with different Poisson’s ratios by assigning different materials. In other words, the patterns that govern the deformation at high temperatures such as the re-entrant cell must be “contained” in the original ground structure. Consequently, our approach is applicable to all ground structures that are statically determinate or overdetermined. This is highlighted by the different choice of ground structures among the three examples in this work, i.e. the 2D square unit cell with a central node, the 3D octet-truss, and the 3D body-centered cubic cell, which are all widely found in literature and are hence chosen as representative cells for this work. Further, the proposed methodology is readily applicable to other types of lattices such as functionally graded structures where the bar thickness varies throughout the structures or different unit cell types, as long as the underlying truss morphology is statically determinate or overdetermined.

Another effect that can be observed in the experiments, especially in the structures “2D Sine” and “3D Sine”, is the buckling of individual bars. As the simulation and the optimization are based on linear-elastic truss elements, local buckling is generally not modelled and must be considered separately, for example by Euler buckling constraints [37]. By comparing the experiments with the simulated structures, it can be seen that only bars that consist of the compliant VW material buckle. At 80 °C, the Young’s modulus of HT is about 36 times higher than the Young’s modulus of VW. Hence, the structural response during the actuation is mostly governed by the stiffer HT members and the influence of the compliant VW members is quite small. The effect of the additional reduction in stiffness of the already compliant members due to the buckling is hence even smaller and can be neglected in the numerical model.

The usage of shape memory polymers allows for deforming the structures temporarily, fixing the deformation, and finally recovering the original shape. As some of the compliant bars buckle during the deformation, these bars are fixed in the buckled shape. Even though the overall topology, i.e. how the bars are connected with each other, does not change, the load-carrying capacity of the structures is influenced by the buckled bars. The buckled bars do not carry as much load as straight bars, which makes the structure in the deformed shape less stiff. However, Fig. 7 shows that the deformed structure is still able to carry a significant amount of load before failure, which can be sufficient depending on the application of the structure. To improve the load-carrying capacity of the deformed structure, a buckling constraint on the normal stress in a bar $\sigma_{{\text{N}}}$ can prevent buckling of the bars during the shape transformation. The critical buckling stress is computed as Euler’s critical stress [37] $\sigma_{{{\text{cr}}}} = (\pi^{2} E_{e} I_{e} /(\left( {0.5l_{e} } \right)^{2} A_{e} )),\quad \forall e = 1 \ldots N_{{\text{b}}}$, with the Young’s modulus $E$, the area moment of inertia $I$, the length of a bar $l$, and the cross-section area of a bar $A$. As a buckling constraint in general reduces the design space, we show the effect of a buckling constraint on the shape morphing problem on the example of the “3D Sine” structure. Figure 8 shows the optimization mean error of the “3D Sine” example with a buckling constraint for different bar radii (blue markers). Each marker represents the mean value of ten optimization runs and the error bars indicate the standard deviation. It can be observed that for larger radii $r \ge 1.75\,{\text{mm}}$, mean errors of about 0.3 mm and below are achieved, which is similar to the values achieved without buckling constraint (black markers). For smaller radii r < 1.75 mm, the mean error increases up to about 2.0 mm. No valid solution can be found for r = 0.6 mm and r = 0.75 mm as the bars are too slender and the prescribed target shape cannot be obtained without buckling. Hence, buckling constraints can help to prevent buckling of members during the shape transformation, which increases the load-carrying capabilities in the deformed state. However, when the bars are too slender, target deformations cannot be achieved at all or only with increasing the shape error.

Figure 8

Optimization mean error of “3D” sine structure with (blue markers) and without (black markers) buckling constraint for different bar radii.

A limitation of the current approach is the heating that is required to enable the deformation of the structures. While heating a full structure evenly is feasible in an experimental setup, for example in an environmental chamber or in a water tank, it is difficult to achieve under real-life conditions. A possible solution to this can be to limit the number of bars that have to be heated to smaller regions, which then can be heated for example by heating wires. This does not require any conceptual changes in the optimization, which makes the whole design framework more versatile.

Another limitation is the assumption of linear elasticity, which limits the magnitude of the achievable target shapes to small deformations. However, it enables the efficient set up of the gradient-based optimization algorithm and makes the method fast and feasible also for large problems with many bars, compared to a nonlinear approach. Hence, the range of deformation is, in the current state of research, a trade-off between efficient modelling and large displacements.

The advantage of our approach compared to other approaches from literature, where individual actuators must be manually programmed for a one-way deformation, lies in the simplified actuation via a single input displacement and the reversibility of the deformation. The optimization setup described in “Methodology” ensures that the deformation is evenly distributed over the whole structure and the original cells remain structurally intact and the efficient optimization problem formulation makes the approach applicable to large-scale structures with at least 3,500 design variables.

Conclusions

In summary, we introduce a novel concept to computationally design lattice structures that can actively change their shape in a reversible way. In contrast to conventional, rigid lattice structures, these structures can temporarily transform from their initial shape to a defined target shape and back. This is achieved by the combination of two different active 3D-printed materials, which have the same stiffness at room temperature but a different stiffness at higher temperatures, and a single, input displacement that controls the global shape of the structures. A gradient-based optimization framework allows for the efficient computational design of large active lattices with 3’500 individual bars. We show that both 2D and 3D structures can be computationally designed, fabricated, and experimentally deformed to achieve the target shapes. Due to the generality of our approach, it can readily be applied to other ground structure types and combinations of boundary conditions and actuation patterns. This makes the approach suitable for many types of applications in aeronautics and automotive, where distinct mechanical properties and shapes for multiple operational conditions are required, for example in active airfoils and shape morphing car panels.

Methodology

Finding a material distribution that yields a specific target shape for prescribed input displacements at 80 °C can be translated to solving a discrete optimization problem. Figure 3A shows the general optimization problem setup based on a 2D lattice structure. A lattice structure in the design domain $\Omega$ is described by the positions of the nodes ${\mathbf{x}}$, how these nodes are connected by $e = 1 \ldots N_{{\text{b}}}$ bars of the length $l_{e}$, the cross-section area $A_{e}$ of each bar, and the material of each bar, represented by the respective Young’s modulus $E_{e}$. To support the structure and to prevent any rigid body modes, a set of Dirichlet boundary conditions is applied. In this example, all nodes on the left side are fixed in the x-direction and in the y-direction, and all nodes on the right side are fixed in the y-direction. A prescribed input displacement ${\mathbf{u}}_{{\text{p}}}$ is applied to all nodes on the boundary $\partial \Omega_{{\text{p}}} \subset \Omega$, which consists of all nodes on the right side in this example. This prescribed input displacement represents the external deformation, for example applied by a pneumatic actuator. The $n_{{\text{T}}}$ target displacements ${\tilde{\mathbf{u}}}_{{\text{T}}}$ represent a sine-wave with an amplitude of $a = 3\,{\text{mm}}$ and are defined on the boundary $\partial \Omega_{{\text{T}}} \subset \Omega$, which includes in total 10 nodes at the top and at the bottom in this example. Since the optimization problem is a problem of finding a discrete material distribution, the design variables are the Young’s moduli $E_{e} \in \left\{ {E_{{{\text{VW}}}} ,E_{{{\text{HT}}}} } \right\}$ of the bars. Each bar can either be made from VW with $E_{{{\text{VW}}}} = 8.5\,{\text{MPa}}$ or from HT with $E_{{{\text{HT}}}} = 305\;{\text{MPa}}$. This model represents the deformation of the structure at a temperature of 80 °C. To avoid only local deformations at large stresses and strains with premature material failure, a constraint on the maximum strain $\varepsilon_{{{\text{max}}}}$ is added. Since the failure stress of HT at 80 °C is much higher than the failure stress of VW, but the failure strain of both materials is approximately $\varepsilon_{{{\text{max}} }} = 0.2$, a strain constraint is more feasible than a stress constraint for this problem. In summary, the optimization problem can be formulated as follows.

$$\mathop {\min }\limits_{E} f\left( E \right) = \mathop \sum \limits_{j = 1}^{{n_{T} }} \left( {u_{{{\text{T}},j}} \left( E \right) - \tilde{u}_{{{\text{T}},j}} } \right)^{2}$$

(1)

$${\text{s.t.}}\quad \varepsilon_{e} \le \varepsilon_{\max } \quad \forall e = 1 \ldots N_{b}$$

(2)

$${\text{with}}\quad E_{e} \in \left\{ {E_{VW} , E_{HT} } \right\}\quad \forall e = 1 \ldots N_{b}$$

(3)

$${\mathbf{Ku}} = {\mathbf{f}}$$

(4)

where ${\mathbf{Ku}} = {\mathbf{f}}$ is the finite element (FE) equation [38] with the stiffness matrix ${\mathbf{K}}$, the displacements ${\mathbf{u}}$, and the loads ${\mathbf{f}}$, and the design variable vector ${\mathbf{E}} = \left[ {E_{1} ,E_{2} , \ldots ,E_{{N_{{\text{b}}} }} } \right]$. For simplicity, all bars are modelled as linear-elastic truss elements. The objective function measures the distance between the target displacements ${\tilde{\mathbf{u}}}_{{\text{T}}}$ and the actual displacements computed by the finite element analysis ${\mathbf{u}}_{{\text{T}}}$, which are evaluated at $n_{{\text{T}}}$ degrees of freedom on the boundary $\partial \Omega_{{\text{T}}}$.

The number of design variables and the number of strain constraints equals the number of bars in the structure, which can quickly become very large for more complex problems or for 3D structures. Previous work shows that solving this type of problem with discrete stochastic solvers such as the Mesh Adaptive Direct Search algorithm (MADS) is only feasible for less than about 200 individual bar [32]. More efficient techniques for solving large-scale optimization problems such as mathematical programming mostly require continuous problem formulations and rely on gradient information. Hence, we relax our problem by making the design variables $E_{e}$ continuous, with the two discrete values of the Young’s modulus of VW and of HT as the lower bound ($E_{{{\text{lb}}}} = E_{{{\text{VW}}}}$) and the upper bound ($E_{{{\text{ub}}}} = E_{{{\text{HT}}}}$), respectively. To guide the optimization towards discrete solutions, a quadratic penalization term, $\varphi$, is added to the objective function:

$$\varphi \left( {\mathbf{E}} \right) = \mathop \sum \limits_{e = 1}^{{N_{b} }} - \frac{{p_{{{\text{max}}}} }}{{\left( {\frac{{E_{{{\text{ub}}}} - E_{{{\text{lb}}}} }}{2}} \right)^{2} }} \times \left( {E_{e} - \frac{{E_{{{\text{ub}}}} - E_{{{\text{lb}}}} }}{2}} \right)^{2} + p_{{{\text{max}}}}$$

(5)

which is zero for $E_{e} = E_{{{\text{lb}}}}$ and $E_{e} = E_{{{\text{ub}}}}$, and has the value $p_{{{\text{max}}}}$ for $E_{e} = \left( {E_{{{\text{lb}}}} + E_{{{\text{ub}}}} } \right)/2$. In combination with the bounds, the penalization term penalizes all values that are not equal to $E_{{{\text{lb}}}}$ or $E_{{{\text{ub}}}}$. Similar to the penalization in material interpolations schemes such as SIMP, a continuation scheme is used to gradually guide the algorithm to “good” solutions and avoid premature local optima [39]. The penalization factor $p_{{{\text{max}}}}$ is zero initially and is increased twice by 0.01, if the objective function value remains below a threshold $f < f_{{\text{p}}} = ({{n_{{\text{T}}} \times a}})({100})$ for 50 iterations. This heuristic adaptation scheme ensures that the penalization is proportional to the magnitude of the target displacements $a$ and the number of specified target displacements $n_{{\text{T}}}$ for each problem, i.e. the problem size and the deformation magnitude. The parameters for the adaptation scheme are obtained heuristically for these examples and may need to be tuned for other problems with different problem sizes and objective functions. The adaptation scheme is well tuned when the penalty value $\varphi$ is close to zero at the end of the optimization and the design variables are close to the bounds, i.e. they are driven to the discrete values.

To avoid costly evaluation of the FE equation and to compute the gradients of the optimization problems efficiently, the FE equation is disaggregated into the static equilibrium equation, the elastic compatibility equation and the constitutive law:

$${\mathbf{Bt}} = {\mathbf{f}}$$

(6)

$${\mathbf{B}}^{{\text{T}}} {\mathbf{u}} = {{\varvec{\epsilon}}}$$

(7)

$${{\varvec{\upsigma}}} = {\mathbf{E\varepsilon }}$$

(8)

with the equilibrium matrix ${\mathbf{B}}$, the member forces ${\mathbf{t}}$, the external loads ${\mathbf{f}}$, the bar elongations ${{\varvec{\epsilon}}}$, the axial strain in the bars ${{\varvec{\upvarepsilon}}}$, and stresses in the bars ${{\varvec{\upsigma}}}$. Note that this formulation is only valid for a truss-based FE discretization. Any point in the solution space that is a solution to these three equations is automatically a valid solution to the original FE equation. With this, the relaxed optimization problem can be reformulated as [40]:

$$\mathop {\min }\limits_{{{\mathbf{E}},{\mathbf{q}}_{{\text{t}}}^{+}, {\mathbf{q}}_{{\text{t}}}^{-}, {\mathbf{u}}}} f_{{\text{r}}} = \mathop \sum \limits_{j = 1}^{{n_{{\text{T}}} }} \left( {{\mathbf{u}}_{{{\text{T}},j}} - {\tilde{\mathbf{u}}}_{{{\text{T}},j}} } \right)^{2} + \varphi \left( {\mathbf{E}} \right)$$

(9)

$${\text{s.t.}}\quad \sum\limits_{(e = 1)}^{{(N_{b} )}} {(q_{{({\text{t}},e)}}^{ + } + q_{{({\text{t}},e)}}^{ - } )} {\mathbf{PC}}_{e} {\mathbf{x}} = {\mathbf{f}}$$

(10)

$$\frac{{A_{e} }}{{l_{e} }}E_{e} {\mathbf{x}}^{{\text{T}}} {\mathbf{C}}_{e} {\mathbf{P}}^{{\text{T}}} {\mathbf{u}} = \left( {q_{{{\text{t}},e}}^{ + } + q_{{{\text{t}},e}}^{ - } } \right)l_{e}^{2}$$

(11)

$$\frac{{\left( {q_{{{\text{t}},e}}^{ + } + q_{{{\text{t}},e}}^{ - } } \right)}}{{E_{e} \varepsilon_{{{\text{max}}}} }} - \frac{{A_{e} }}{{l_{e} }} \le 0$$

(12)

$${\text{with}} \quad E_{e} \in \left[E_{\text{VW}}, E_{\text{HT}}\right]$$

(13)

$\forall e = 1 \ldots N_{{\text{b}}}$. In this formulation, Eqs. (6)–(8) are combined into Eqs. (10) and (11). Equation (12) represents the strain constraint. The non-negative, normalized member forces $q_{{{\text{t}},e}}^{ + } + q_{{{\text{t}},e}}^{ - } = q_{{{\text{t}},e}} :\, = t_{e} /l_{e}$ and the displacements ${\mathbf{u}}$ are directly treated as design variables, which avoids the costly inversion of the stiffness matrix ${\mathbf{K}}$. The matrices ${\mathbf{P}}$ and ${\mathbf{C}}_{e}$ are used for an efficient implementation of the equilibrium matrix ${\mathbf{B}}$. A detailed description of the full optimization problem as well as the gradients with respect to the design variables can be found in the Supplementary Information (SI). The optimization problem is solved using IPOPT, a state-of-the-art interior point line search filter method for solving large-scale, nonlinear, constrained, bounded optimization problems [41] and the MATLAB interface of the OPTI Toolbox [42].

Even though the relaxed optimization model penalizes intermediate design variable values, i.e. intermediate Young’s moduli, it does not ensure that all design variables are discrete at the end of the optimization, which can be the case when the solver converges to a local minimum. Hence, to ensure a discrete material distribution that can be fabricated using the chosen 3D printing material set, a “postprocessing” step is required. One possibility for the postprocessing step can be to apply a projection scheme where, for example, all values closer to the lower bound or the upper bound are set equal to the respective bounds. If the penalty term in the initial relaxed optimization problem is tuned well as described previously, only few design variables are affected by the projection. However, especially for small problems with fewer design variables the quality of the solution can be affected significantly. A possibility to simultaneously ensure a discrete material distribution and consider the quality of the solution is to use a stochastic optimization method as a postprocessing step, with the solution from the initial IPOPT optimization as a starting point. Since the stochastic optimization starts from an “almost discrete” solution, the complexity of the problem is already significantly reduced compared to a full discrete optimization, as for example described in literature [32]. This makes it numerically cheaper than a full discrete optimization, but more accurate than a simple projection. Here, the MATLAB Genetic Algorithm (GA) is chosen for the postprocessing step since it can handle discrete, large-scale optimization problems. Hence, in the postprocessing step, the discrete problem described by Eqs. (1)–(4) is solved directly starting from the solution from the relaxed optimization problem described by Eqs. (9)–(13). All bar materials remain design variables but can now only take on discrete values and all intermediate values are projected to the closest discrete value. While this introduces more randomness to the optimization, the general quality of the solution is preserved as the solution from the first step with a low objective function is used as a starting point. Hence, the random variations in the postprocessing step help to identify similar configurations that might have better objective function values than the projected solution after the first step. The finite element equation in Eq. (4) is directly solved by an inhouse truss FEM implementation written in MATLAB. The standard MATLAB GA options are used, except for the maximum number of generations and the population size, which are both set to 150.

In summary, every optimization run consists of two steps: (1) solving the relaxed optimization problem with IPOTP, (2) postprocessing the solution from step 1 by solving the discrete optimization problem with a GA. Since the optimization problem is nonlinear and non-convex, the solution of step 1 depends on the starting point and finding a global optimum is not guaranteed. The GA in step 2 further introduces a random component to the problem. To estimate the reliability of our optimization approach, each example problem is solved 10 times and the arithmetic mean and the standard deviation for both steps are reported in Table 1. Step 1 of the optimization framework is solved with random starting points such that all design variables lie within the bounds. For all three example structures, the best solution after both steps is selected out of the 10 optimization runs, fabricated and experimentally tested.