Computational design and fabrication of active 3D-printed multi-state structures for shape morphing

This work introduces a new form of multi-state fabrication to encode multiple motions into a single actuated object by exploiting the heat sensitivity of shape memory polymers (SMPs). Popular approaches to 4D printing exploit the special molecular arrangement of SMPs, which allow deformations in the rubbery state to be fixed by cooling the material below the glass transition temperature ${T_{\text{g}}}$ while holding it in its deformed shape [10]. This is often referred to as the programming step in literature [11]. By re-heating the SMP above its ${T_{\text{g}}}$, the undeformed initial state can be recovered. This is normally achieved via heating in a water bath or with warm air. Recovery to the original state is driven by low magnitude internal forces and repeated actuations or actuations to different states require re-execution of the tedious programming step. Most 4D printed structures today leverage the active properties of their base materials to change their shape [5]. The shape-memory effect, enabled by controlled temperature changes, and simple affine transformations such as shrinking or expansion are translated to more complex shape changes by means of topological and geometric design [12]. This involves 1D rods that deform into 2D objects [13], planar 2D structures that deform into complex 3D structures [14], and 1D rods that assume 3D shapes [15]. The combination of multiple materials with multiple or reversible actuation modes further enables advanced deformations of structures. Composites of multiple active materials or of active and passive materials, such as other polymers, evoke local stress and strain anisotropies and can trigger global deformations [11, 16]. Similar behavior can be observed by printing structures with pre-strains, which can be released by adding a temperature stimulus and in turn trigger deformation [17]. Other studies use numerical simulations and analytical models to investigate fundamental deformation modes of active joints, bars, and hinges such as bending and curling [6, 15, 18–20]. However, these structures are often limited to a single, pre-programmed shape, require global heating for example in a water bath, or can only be actuated once without re-programming.

The approach presented in this work significantly differs from many standard 4D printing approaches as the motion of the structures is not generated by the recovery forces of the SMP, but by an independent input actuation. Rather than an ambient heat source, copper wires are used to locally control the heating of the structures. This creates local material stiffness gradients, which alter the global deformation mode for a given driving actuation. By selectively heating different parts of the structure, different end states can be generated despite starting from the same initial state of a single printed structure. This approach confers two advantages. First, the structures can transition between states relatively quickly as the transition speed is determined by the speed of the applied driving actuation. And second, structures can repeatedly transition between states without re-programming as the state transitions are encoded once at the moment of fabrication and remain persistent thereafter.

The key computational problem in the multi-state fabrication design is one of material placement and heating: how can material be arranged and heated inside the design domain such that objects can achieve the correct deformed shapes? Hence, the input to the optimization framework is a design domain ${{\Omega }}$ that describes the initial shape of the multi-state object, a given base material (including mechanical and thermal properties), a set of boundary conditions, a single driving actuation, and a set of target deformations. The output is a material distribution and a heating pattern, which enable the desired target deformations for the prescribed input actuation and boundary conditions. To obtain structures that can be fabricated and actuated, several requirements are imposed on the optimization scheme. As such, the method must:

• match multiple target states under a given set of boundary conditions.
• Have a binary material distribution.
• Form a single, connected region.
• Limit material feature size to prevent regions that are too small to be printed or actuated.

The overall design approach comprises several steps and is illustrated in the following paragraphs using the example of a 'Gingerbread Man' that can lower its left arm or lower its right arm by pushing its head, schematically shown in figure 1.

Zoom In Zoom Out Reset image size

To begin, a standard single-state deformation design problem is considered. The goal is to minimize the pose error $J$, which measures the distance between the displacements $u$ of a body under prescribed boundary conditions and the target displacements ${u_{\text{T}}}$ evaluated on the target boundary ${{{\Gamma }}_{\text{T}}}$, with the material density distribution ${\boldsymbol{{\rho }}}$ of the body as the design variable. The standard ${L_2}$-norm is used to measure the error between $u$ and ${u_{\text{T}}}$. This yields the following design optimization problem with the objective function $J$:

$$\begin{align}{{\text{min}} _{\boldsymbol{\unicode{x03C1} }}}J\left( {\boldsymbol{\unicode{x03C1} }} \right) = \left\| {{u_{\text{T}}} - u\left( {\boldsymbol{\unicode{x03C1} }} \right)} \right\|.\end{align} \tag{ 1 }$$

The structure is fixed at the boundary ${{{\Gamma }}_{\text{f}}}$. This basic design optimization problem essentially describes the design of a compliant mechanism as it is often found in literature [21], where an input displacement or force is transferred to an output deformation or force by means of elastic deformation of the structure.

This formulation is easily extended to include multiple deformed states by considering the summation of individual pose errors, ${J_j}$, one for each set of $k$ target deformations:

$$\begin{align}\begin{gathered} {\text{min}}{_{\boldsymbol{\unicode{x03C1} }}}\sum\limits_{j = 1}^k {{J_j}\left( {\boldsymbol{\unicode{x03C1} }} \right)} \\ {\text{with}}{\,J_j}\left( {\boldsymbol{\unicode{x03C1} }} \right) = \left\| {{u_{{\text{T}},j}} - {u_j}\left( {\boldsymbol{\unicode{x03C1} }} \right)} \right\| \\ \end{gathered}.\end{align} \tag{ 2 }$$

Next, it must be ensured that the deformations of the structure are physically valid. This is achieved by adding a constraint to enforce static equilibrium of the structure. Assuming a linearly elastic constitutive model, this constraint can be discretized using the finite element method. For a single target state, this yields:

$$\begin{align} {\text{min}}{_{\boldsymbol{\unicode{x03C1} }}}J\left( {\boldsymbol{\unicode{x03C1} }} \right)& = \left\| {{{\mathbf{L}}_{\text{u}}}\left( {{u_{\text{T}}} - u\left( {\boldsymbol{\unicode{x03C1} }} \right)} \right)} \right\| \nonumber \\ {\text{for}}\;{{\mathbf{K}}_{\boldsymbol{\unicode{x03C1} }}}{\mathbf{u}} &= {\mathbf{f}} \end{align} \tag{ 3 }$$

where ${{\mathbf{K}}_{\boldsymbol{\rho }}}$ is the stiffness operator. The binary diagonal matrix ${{\mathbf{L}}_{\text{u}}}$ restricts $J$ to just the degrees of freedom that correspond to the specified target deformations on the boundary ${{{\Gamma }}_{\text{T}}}$. Hence, the diagonal entries of $J$ are simply the indicator function for ${{{\Gamma }}_{\text{T}}}$. The constraint ensures that the deformed state is in static equilibrium.

Optimizing only the material densities ${\boldsymbol{\unicode{x03C1} }}$ is not sufficient to enable multi-state deformation resulting from a single input actuation. Rather the optimization must include the possibility to take advantage of the temperature varying stiffness of the SMP material. This is accomplished by introducing an additional optimization variable, $\eta \,:{{\Omega }} \to \left[ {0,1} \right]$, which encodes the thermal state of the material at each point in the object. Here, $\eta = 0$ represents the material at room temperature, while $\eta = 1$ indicates that the material is heated above the ${T_{\text{g}}}$. In combination with the material density, ${\boldsymbol{\unicode{x03C1} }}$, which controls the presence or absence of material [22], the following material interpolation scheme is used for computing the Young's modulus:

$$\begin{equation}E\left( {\rho ,\eta } \right) = \mathop \int \limits_{{\Omega }}^{} \rho \left[ {\left( {1 - \eta } \right){E_{{\text{min}}}} + \eta {E_{{\text{max}}}}} \right]\end{equation} \tag{ 4 }$$

where ${E_{{\text{min}}}}$ and ${E_{{\text{max}}}}$ are the minimum and maximum stiffness of the material, which correspond to temperatures above and below the ${T_{\text{g}}}$, respectively. This describes the design optimization problem as a combined topology and material optimization formulation, where every element has the two design variables $\rho$ and $\eta$. To prevent the stiffness matrix to become singular, the minimum density value is set to $\rho = 0.001$ as commonly described in literature [22].

To ensure that the optimization algorithm generates structurally sound and connected topologies, a compliance term $C\left( {{\boldsymbol{\unicode{x03C1} }},{\boldsymbol{\eta }}} \right)$ is introduced to the objective function. This approach is similar to others found in literature [21, 23, 24] but adapted to the particular displacement-driven use case in this work. This additional term is required as many configurations exist that are kinematically valid and produce the correct target deformations, but are statically not necessarily feasible, as they for example include very thin connections or unconnected regions. To bias the optimization towards solutions that fulfill both criteria, the observation that connected, structurally sound topologies are in general stiffer than unconnected ones can be leveraged. Hence, additional compliance terms are introduced, one for each target deformation state, that measure the compliance in response to unit reaction forces generated by ${u_{\text{T}}}$. These reaction forces ${f_{\text{C}}}$ act in the opposite direction of ${u_{\text{T}}}$ and have normal magnitude. These additional terms maximize the stiffness between the target boundary ${{{\Gamma }}_{\text{T}}} = {{{\Gamma }}_{\text{C}}}$ and the fixed boundary ${{{\Gamma }}_{\text{f}}}$, ensuring that both regions are connected, and no unconnected regions are evolving. To comply with the overall minimization problem formulation, a minimum compliance formulation is used, which is equivalent to maximizing the stiffness. The resulting compliance energy $C$ is obtained by solving another linear elasticity problem for each target deformation state to compute ${u_{\text{c}}}$:

$$\begin{align} C\left( u \right) &= \mathop {{{\mathop \int \nolimits }}}\limits_{{{{\Gamma }}_{\text{c}}}}^{} {f_{\text{c}}}{u_{\text{c}}}\\ {\text{with }}{K_{\rho ,\eta }}{u_{\text{c}}} &= {f_{\text{c}}} . \nonumber \end{align} \tag{ 5 }$$

Finally, it must be ensured that the output of the design optimization algorithm is a binary set of material/empty space assignment, i.e. $\rho \in \left\{ {0,1} \right\}$. This is achieved by three modifications of the original optimization problem. First, the SIMP penalization is used with the parameter $p = 3$ to penalize intermediate densities and thermal states [25]. This amounts to a simple modification of equation (4) to

$$\begin{equation}E\left( {\rho ,\eta } \right) = \mathop {{{\mathop \int \nolimits }}}\limits_{{\Omega }}^{} {\rho ^{\,p}}\left[ {\left( {1 - {\eta ^{\,p}}} \right){E_{{\text{min}}}} + {\eta ^{\,p}}{E_{{\text{max}}}}} \right].\end{equation} \tag{ 6 }$$

However, this adaptation is not enough to promote binary designs, as it generally promotes 'efficient' material usage for compliance-based optimization problems. As minimizing a compliance term is only a part of the problem formulation here, additional adjustments are necessary to promote fully solid-void-designs. Hence, a regularization term $R\left( \xi \right)$ is additionally introduced that penalizes the intermediate design variables $\xi \in \,\left] {{\xi _{{\text{min}}}},{\xi _{{\text{max}}}}} \right[$. The regularization takes the form of a quadratic function:

$$\begin{equation}R\left( \xi \right) = - \frac{{{R_{{\text{max}}}}}}{{{{\left( {\frac{{\left( {{\xi _{{\text{max}}}} - {\xi _{{\text{min}}}}} \right)}}{2}} \right)}^2}}}{\left( {\xi - \frac{{{\xi _{{\text{max}}}} + {\xi _{{\text{min}}}}}}{2}} \right)^2} + {R_{{\text{max}}}}\end{equation} \tag{ 7 }$$

since the regularization term is $0$ only at the bounds ${\xi _{{\text{min}}}}$ and ${\xi _{{\text{max}}}}$, a penalty is added to all intermediate values. If any intermediate values remain after the regularization, a projection $P\left( \xi \right)$ is applied. For some scalar value $\xi$, the projection maps all intermediate values $\xi \in \,\left] {{\xi _{{\text{min}}}},{\xi _{{\text{max}}}}} \right[$ to the space of binary values $\left\{ {{\xi _{{\text{min}}}},{\xi _{{\text{max}}}}} \right\}$

$$\begin{equation}P\left( \xi \right) = \left\{ {\begin{array}{*{20}{c}} {{\xi _{{\text{min}}}}\,\,\,\,\,\,if\,\,\xi \leqslant {\xi _{\text{t}}}} \\ {{\xi _{{\text{max}}}}\,\,\,otherwise} \end{array}} \right.\end{equation} \tag{ 8 }$$

here, the threshold ${\xi _{\text{t}}} \in \left( {{\xi _{{\text{min}}}},{\xi _{{\text{max}}}}} \right)$ is chosen such that the target deformation objective $J$ is minimized. This yields the optimization problem formulation that describes the computational design problem for $k$ target deformation states:

$$\begin{align}\begin{array}{l} \mathop {\min }\limits_{\rho ,{\mathbf{\eta }}} {\mkern 1mu} \sum\nolimits_{j = 1}^k {{J_j}\left( {\rho ,{{\mathbf{\eta }}_j}} \right)} + \alpha \sum\nolimits_{j = 1}^k {{C_j}\left( {\rho ,{{\mathbf{\eta }}_j}} \right)} \\ + \sum\nolimits_{j = 1}^k {\sum\nolimits_{i = 1}^n {{R_i}\left( {\rho ,{{\mathbf{\eta }}_j}} \right)} } \\ {\textrm{ s}}.{\textrm{t}}.{\textrm{ }}{{\mathbf{K}}_{\rho ,{\eta _j}}}{{\mathbf{u}}_j} = {{\mathbf{f}}_j},{\textrm{ }}\forall {\textrm{j}} = 1 \ldots k\;\\ {{\mathbf{K}}_{\rho ,{\eta _j}}}{{\mathbf{u}}_{{\textrm{c}},j}} = {{\mathbf{f}}_{{\textrm{c}},j}},{\textrm{ }}\forall {\textrm{j}} = 1 \ldots k \end{array} \end{align} \tag{ 9 }$$

where $\alpha$ is a weighting parameter as commonly used in multi-objective optimization problems, which is set to $\alpha = 1$ in this work.

This section describes the finite element discretization, which lies at the heart of the proposed method. Although the method relies on standard linear elasticity, and displacement and traction boundary conditions are applied in the normal way, the geometric discretization is unique in that it is based on adaptive, volume constrained power diagrams. Typical approaches to topology optimization rely on extremely refined meshes to capture geometric details, which can lead to a variety of artifacts [22]. This specifically includes checkerboarding, hinging, and very small features that negatively influence the fabricability [26]. To avoid these complications and to facilitate the actuation of the structures in this work, this approach relies on constructing a comparably small number of polygonal elements that can adapt themselves to find an optimized arrangement. Every polygon represents one individual cell with a dedicated material density and individual heating scheme. Hence, the size of the polygons is directly related to the size of the structural features within the design domain. The advantages of the polygon-based approach are highlighted in figure 2. The structure on the left is discretized by individual triangular elements, where every finite element has its own density and heating pattern. On the right, a structure partitioned by polygonal cells is shown, where the density and heating pattern is defined for every cell. While both structures achieve the same deformation behavior, the polygonal meshing reduces the overall number of variables, explicitly avoids high-resolution artifacts such as the checkerboarding, and hence improves the fabricability of the structure. Finally, making the geometry of the cells adaptive and including them as design variables in the optimization routine allows the underlying mesh to change its geometry and tune the shape of the resulting structures independent of the density distribution.

Zoom In Zoom Out Reset image size

In this work, power cells are used to generate the adaptive polygonal partitioning of the design domain ${{\Omega }}$. Power cells are a generalized version of Voronoi cells. Voronoi cells in 2D describe convex polygons that are based on a set of points in the Euclidean plane. In contrast to Voronoi cells, power cells can change their geometry and size by adding specific weights to the cells. Figure 3 shows a set of power cells and the corresponding nomenclature. A power diagram [9] partitions a domain ${{\Omega }}$ into cells, where every power cell ${\mathcal{V}_i}$ is defined by the respective power diagram site ${x_i}$ and weight ${w_i}$ via

$$\begin{align}\mathcal{V}_i^w = \left\{ {x \in {{\Omega }} |\ \ \|{ }x - {x_i}\|^2 - {w_i} \leqslant \|\,x - {x_j}\|^2 - {w_j},\,\forall j} \right\}\end{align} \tag{ 10 }$$

Zoom In Zoom Out Reset image size

The two cells ${\mathcal{V}_i}$ and ${\mathcal{V}_j}$ are neighbors if the intersection ${\mathcal{V}_i}\mathop \cap \nolimits^ {\mathcal{V}_j}$ is not empty (an edge in 2D or a planar polygon in 3D). The boundary between the two cells is denoted as $e_{ij}^*$. The Euclidean distance between ${x_i}$ and ${x_j}$ is $\left| {{e_{ij}}} \right|$ and the area of a power cell is ${V_i}$. Note that a power diagram with equal weights of all cells equals the Voronoi diagram with the same sites. The well known dualism between Voronoi and Delaunay cells also holds for power diagrams, where the dual cells constructed from connections between power sites orthogonal to the cell boundaries are called weighted Delaunay cells [27].

To increase the generality of the proposed method and include the design of structures with non-convex boundary shapes, the construction of power diagrams on non-convex domains must be adapted as power cells by definition are convex. For a given non-convex domain ${{\Omega }}$, this is done by first constructing the normal power diagram with $n$ randomly initiated sites inside of the domain. Next, all cells that intersect with ${{\Omega }}$ are cut by applying the Sutherland–Hodgeman algorithm [28]. For non-convex domains, this can result in the formation of non-convex power cells, from which all relevant parameters such as ${V_i}$, $e_{ij}^*$, and ${e_{ij}}$ can then be computed.

For structural computations, a triangular finite element mesh is then constructed from the clipped cells by connecting the two ends of every edge of a power cell ${\mathcal{V}_i}$ with the power cell site ${x_i}$. The cell boundary vertices and the cell sites form the vertices $\chi$ of the finite element mesh. For non-convex cells, not all boundary vertices can be connected to the cell sites. To build a triangular mesh in these cells without introducing additional points and edges, a classic Delauney triangulation [29] is employed that only uses cut-polygon edges and the cell site. The dotted red lines in figure 3 show the triangular mesh constructed from convex and clipped non-convex power diagram cells. In the optimization framework, the same material density and heating scheme is assigned to all triangular finite elements in one cell. Hence, the number of cells directly links to the number of design variables.

This meshing approach further simplifies the calculation of derivatives compared to other meshing strategies since the finite element mesh vertices are now implicitly linked to the construction of the (non-)convex power cells.

A key advantage of power diagrams over Voronoi diagrams is that the cell sizes can be changed without having to move the cell sites. This is useful as by constraining the size of the individual cells, a minimum feature size can be enforced. This is implemented by adjusting the weights ${w_i}$ of the power cells. It was shown by [30] that finding a set of weights that enforces the cell volumes ${V_i}$ to match a given set of target cell volumes ${V_{{\text{t}},i}}$ is equivalent to minimizing the smooth energy

$$\begin{equation}\mathcal{E}\left( {X,W} \right) = {\sum\nolimits_i {\mathop \int \limits_{{\mathcal{V}_i}}^{}} {\left\| {x - {x_i}} \right\|} ^2}dx - \sum\nolimits_i {{w_i}\left( {{V_i} - {V_{{\text{t}},i}}} \right)}. \end{equation} \tag{ 11 }$$

De Goes et al [31] show that minimizing this energy is equivalent to solving an optimal transport problem, a well-studied problem in computer graphics, with volume constraints. The set of weights $W$ is found by updating the power diagram via Newton iterations and the derivatives have simple, closed-form expressions. It is further shown that this directly relates to finding a centroidal power diagram. Centroidal power diagrams are power diagram where each site ${x_i}$ coincides with the centroid, i.e. the center of mass, of each cell.

Since a volume-constrained power diagram is unique for a given set of target volumes ${V_{{\text{t}},i}}$, the geometry of the power cells can be inversely controlled by prescribing different target volumes. This behavior can be used in the overall optimization framework proposed here by treating the target volumes as design variables to allow the power diagram, and hence the finite element mesh, to adapt during the optimization. Limiting the minimum and maximum cell volume to ${V_{{\text{min}}}} \leqslant {V_i} \leqslant {V_{{\text{max}}}}$ gives direct control over the feature size of the resulting geometry, enabling smooth, fabrication-ready structures. The algorithm to produce centroidal, volume-constrained power diagrams based on [31] is briefly summarized in algorithm 1. The threshold for the gradient norm, which directly measures how close the sites are to the actual centroid of the polygons, is computed via ${\in _{{\text{VCPD}}}} = {10^{ - 4}}{L_{{\text{domain}}}}{V_{{\text{mean}}}}\sqrt {8n}$. This ensures that the mean error between the cell sites and the cell centroids is related to the overall size of the problem and is on average not larger than ${10^{ - 4}}{L_{{\text{domain}}}}$, where ${L_{{\text{domain}}}}$ is the maximum length of the domain ${{\Omega }}$. To ensure that the individual cell volumes ${V_i}$ always add up to the total domain volume ${V_{{\text{total}}}}$, the relative cell volume $\phi$ is introduced. For a given set of variables ${\phi _i}$, the physical cell volumes ${V_i}$ can be computed via

$$\begin{align}{V_i} = \frac{{{\phi _i}}}{{{{\mathop \sum \limits }}_{i = 1}^n{\phi _i}}}{V_{{\text{total}}}}.\end{align} \tag{ 12 }$$

. 

Algorithm 1. Centroidal volume-constrained power diagram.
Data: domain ${{\Omega }}$, target volumes ${V_{\text{t}}}$, number of cells $n$, initial positions of sites $X$
Result: Centroidal volume-constrained power diagram
1 while  ${\|\nabla _X}\varepsilon\| \geqslant {\in _{{\text{VCPD}}}}$  do
2  Enforce volume constraint;
3  Lloyd-step or gradient-descent step;
4 Update power diagram;
5 end

The power diagrams are integrated into the topology optimization by explicitly updating the power diagram at the beginning of every optimization step according to the relative cell volume design variables $\phi$. After that, the finite element mesh is constructed and the objective function is computed. The optimization problem described in equation (9) is modified to become

$$\begin{align} & \mathop{{\text{min}}}\limits_{{\boldsymbol{\unicode{x03C1} }},{\boldsymbol{\eta }},\boldsymbol{\phi} }\,\mathcal{F} = \sum\limits_{j = 1}^k {{J_j}\left( {{\boldsymbol{\unicode{x03C1} }},{{\boldsymbol{\eta }}_j},\boldsymbol{\phi} } \right)} + \alpha \sum\limits_{j = 1}^k {{C_j}\left( {{\boldsymbol{\unicode{x03C1} }},{{\boldsymbol{\eta }}_j},\boldsymbol{\phi} } \right)}\nonumber \\ & \qquad + \sum\limits_{j = 1}^k {\sum\limits_{i = 1}^n {{R_i}\left( {{\boldsymbol{\unicode{x03C1} }},{{\boldsymbol{\eta }}_j}} \right)} } \nonumber\\ & {\text{s}}.{\text{t}}.\ \ \ {{\mathbf{K}}_{\rho ,{\eta _j},\phi }}{{\mathbf{u}}_j} = {{\mathbf{f}}_j},{\text{ }}\forall {\text{j}} = 1 \ldots k \nonumber\\ & {{\mathbf{K}}_{\rho ,{\eta _j},\phi }}{{\mathbf{u}}_{{\text{c}},j}} = {{\mathbf{f}}_{{\text{c}},j}},{\text{ }}\forall {\text{j}} = 1 \ldots k \nonumber\\ & \mathop \sum \limits_{i = 1}^n {V_i}\left( \boldsymbol{\phi} \right) = {V_{{\text{total}}}} \nonumber \\ & {\text{with}}\,\,\,{J_j} = \left\| {{{\mathbf{L}}_{u,j}}\left( {{u_{{\text{T}},j}} - {u_j}} \right)} \right\|,\,\,\,\forall j = 1 \ldots k \nonumber \\ & \qquad {C_j} = f_{{\text{c}},j}^T{u_{{\text{c}},j}},\,\,\,\,\,\forall j = 1 \ldots k\nonumber \\ & \qquad {R_i}\left( {{\boldsymbol{\unicode{x03C1} }},{{\boldsymbol{\eta }}_j}} \right) = {R_i}\left( {\boldsymbol{\unicode{x03C1} }} \right) + {R_i}\left( {{\eta _j}} \right),\,\,\,\,\,\,\forall j = 1 \ldots k,\,\forall i = 1 \ldots n \end{align} \tag{ 13 }$$

to minimize this objective, the gradients of the power diagram mesh with respect to the relative power cell volumes $\boldsymbol{\phi}$ are required, which are provided in the supplementary information. The full optimization algorithm is summarized in algorithm 2.

. 

Algorithm 2. Multi-state topology optimization.
Data: domain ${{\Omega }}$, boundary conditions, $k$ target deformations ${u_{{\text{T}},j}}$, number of power diagram cells $n$
Result: binary material density distribution; heating scheme to achieve different target deformations with asingle structure
1 initialize power diagram sites $X$ with $n$ random points inside the domain ${{\Omega }}$;
2 initialize normalized cell sizes with ${\phi _i} = 1$;
3 initialize cell densities with ${\rho _i} = 0.5$;
4 initialize temperature with ${\eta _i} = 0.5$;
5 initialize ${R_{{\text{max}}}} = 0$;
6 while  ${\|\nabla _{\rho ,\eta ,\phi }}{\mathcal{F}\|_\infty } \geqslant {10^{ - 6}}$ and ${n_{{\text{iter}}}} \leqslant {n_{{\text{iter}},{\text{max}}}}$  do
7   compute volume-constrained power diagram;
8   compute derivatives;
9 update design variables;
10 end
11 increase Rmax
12 while  ${\|\nabla _{\rho ,\eta ,\phi }}{\mathcal{F}\|_\infty } \geqslant {10^{ - 6}}$ and ${n_{{\text{iter}}}} \leqslant {n_{{\text{iter}},{\text{max}}}}$  do
13  compute derivatives;
14 update design variables;
15 end
16 apply projection scheme to obtain binary values;

The two-step approach allows the algorithm to converge to a general solution first without applying the regularization. Such continuation schemes where intermediate densities are not fully penalized right from the beginning are common in topology optimization [22]. Hence, the first step in the optimization is used to find a find a structural configuration of power cells along with a material distribution that minimizes the first two terms of the objective function. In a second step, the regularization term is added, but the power diagram is fixed. This step penalizes intermediate material densities and intermediate temperatures and promotes discrete solutions. MATLAB's fmincon interior point algorithm with a limited-memory broyden–fletcher–goldfarb–shanno (L-BFGS) Hessian approximation is used in both steps to solve the optimization problem. If any non-discrete material densities or temperature values remain, the projection scheme is applied. We use a standard, linear triangular element approach often referred to as constant strain triangle elements or T3 [32], to simulate the deformations of our structures. The correctness of our implementation was validated by comparisons with the commercial FE software Abaqus 6.14 (Dassault Systèmes Simulia Corp., Providence, RI, USA).

To fabricate the multi-state structures resulting from the optimization framework, a Stratasys Objet500 Connex3 3D printer (Stratasys, Eden Prairie, MN, USA) is used. Multi-state fabrication can be achieved by exploiting the differences in Young's modulus $E$ and glass transition temperature ${T_{\text{g}}}$ between the most rigid material VeroWhite+ (VW) and the most compliant material Agilus30 (AG). VW has a glass transition temperature of about ${T_{\text{g}}} \approx 65^{\,\circ} {\text{C}}$ and hence is rigid at room temperature ($23^{\,\circ} {\text{C}}$) with $E \approx 2100{\text{ MPa}}$. Upon heating above its ${T_{\text{g}}}$ to $70^{\,\circ} {\text{C}}$, the Young's modulus is reduced to about $8{\text{ MPa}}$. AG has a glass transition temperature of ${T_{\text{g}}} \approx 15^{\,\circ} {\text{C}}$ and is at room temperature already in its rubbery state with $E \approx 0.8{\text{ MPa}}$. Upon heating it to $70^{\,\circ} {\text{C}}$, the Young's modulus is further reduced to about $0.2{\text{ MPa}}$. The multi-state approach requires a material that has a large change in stiffness at elevated temperatures but is compliant enough so that it can be deformed by hand at room temperature. As this is not directly achieved by any of the two base materials or any of the off-the-shelf intermediate, digital materials provided by the 3D printer, a manually programmed voxel-based dithering algorithm is used in this work. By dithering the two base materials at a voxel level [33, 34], the mechanical and thermal properties of the resulting material can be tuned.

To produce voxelized materials, the solid regions of the printable objects are divided into cuboids with edge lengths of $0.25 \times 0.25 \times 0.40{\text{ m}}{{\text{m}}^{\text{3}}}$. Each cuboid represents one voxel and is randomly assigned one of the two base materials (VW or AG), where the overall percentage of both materials can be prescribed. The material properties of the resulting dithered materials are characterized in tensile tests according to the ASTM norm D638-14 [35]. Figures 4(A) and (B) show the stress–strain curves of different dithered materials with 10%, 30%, 50%, 70%, and 90% AG at room temperature (A) and at $70^{\,\circ} {\text{C}}$ (B), respectively. With increasing percentage of AG, the dithered material becomes significantly more compliant in both states. Figure 4(C) shows the Young's modulus at both temperatures for the different dithering ratios. Since the Young's modulus at high temperatures (red dotted line) is comparably low and does not change significantly over the temperature range as both materials are in their rubbery state, the absolute stiffness reduction by heating of a dithered material can be directly controlled by changing the mixture ratio. For the multi-state structures in this work, a mixture ratio of 50% is chosen, where the Young's modulus at room temperature of ${E_{{\text{RT}}}} \approx 120{\text{ MPa}}$ drops to ${E_{70}} \approx 2.9{\text{ MPa}}$ at $70^{\,\circ} {\text{C}}$. As we only consider discrete initial and target states, the thermo-mechanical behavior of the material and structures during the heating process does not have to be considered in order to design structures with target deformation modes. This simplifies our approach compared to other methods where extensive thermo-mechanical modeling is required to describe the structural behavior at all times.

Zoom In Zoom Out Reset image size

Figure 4(D) schematically shows the experimental set up for the 'Gingerbread Man' structure. Once printed, loops of coated copper wire with a diameter of $0.15{\text{ mm}}$ are used to differentially heat the structures. For every target deformation state, one continuous wire is used. The wire loops are fixed with super glue to the structure and are heated to about $70^\circ {\text{C}}$ by applying a current of about $0.9{\text{ A}}$. For actuation, the structures are fixed to a base plate with metal nails and out of plane deformations are prevented by an acrylic plate on top of the structures. The input actuation displacement is manually applied with an acrylic bar. To produce the results predicted by the simulation and make them comparable, we fix the input actuation displacement at the same magnitude as in the simulation and compare the resulting, discrete states of the structures.