Level Set Topology Optimization of Printed Active Composites

The emerging paradigm of PACs builds upon the capability of advanced 3D printers to place multiple materials precisely in a 3D design space [1,2]. PACs derive their behavior and functional performance from the heterogeneous distribution of multiple photocured polymers where some or all of them may exhibit an active response, such as a thermally induced SM behavior. While the potential of PACs for self-actuating and self-assembling structures was recently demonstrated in Ref. [3], only simple PAC designs with conventional fiber architectures were considered. This study aims at exploiting the high, micrometer resolution of multimaterial printers, which allow the design of PACs with complex material layouts, not limited to conventional fibers–matrix arrangements. We will demonstrate that optimally arranging active and passive phases enables PACs to assume desired complex shapes upon temperature change.

The design of PACs is a nontrivial inverse problem of finding the spatial arrangement of active and passive material phases. The active material undergoes complex thermomechanical transformations during the SM cycle, which includes the change of storage modulus and partial stress relaxation upon phase change. The final shape of the PAC is a result of the SM cycle, the mechanical interplay between active and passive phases, and the application of programing loads [4–6]. To account for these phenomena via appropriate mechanical models, we introduce a computational design optimization approach that finds the material layout such that the activated PAC assumes a target shape. To search for the optimal design among a large set of possible designs, including nonintuitive arrangements of materials, we adopt a topology optimization approach.

The two main approaches for structural topology optimization are density methods and LSMs. Adopting a density method for a two-phase material layout problems, a fictitious porous material is introduced and the density is considered a design variable which can continuously vary between ρ = 0, representing phase “1”; and ρ = 1, representing phase “2.” The material properties are interpolated as functions of the density, using some form of homogenization method or via analytical models. The latter approach is used, for example, in the solid isotropic material with penalization (SIMP) method [7,8]. SIMP and other density methods have been successfully applied to a broad range of engineering design problems [9,10]. For problems where the material phases undergo a complex constitutive behavior, as it is the case for PACs, density methods face the challenge of finding appropriate material interpolation rules and penalization schemes that lead to a “0–1” material distribution. The presence of fictitious material with intermediate densities does not only affect the accuracy in predicting the performance of the fabricated PAC but also requires additional postprocessing steps for extracting the geometry of the phase boundaries. The latter step is required to provide 3D printers with a crisp geometry description and may add additional errors in the design process. To overcome these issues, we adopt a level set optimization approach for designing PACs.

LSMs provide a crisp description of the phase boundaries which are defined by an isocontour of a level set function (LSF) [11]. Smooth variations of the LSF lead to shape changes and may also alter the topology of the design. The explicit definition of the phase boundaries facilitates the fabrication of the design, without introducing errors when exporting the geometry from the computational model to a fabrication tool, such as a 3D printer. To preserve the crispness of the geometry model for predicting the response of the PAC by finite elements, we combine the LSM with the XFEM [12]. In contrast to ersatz material approaches, which use material interpolation for mapping the LSF onto a finite-element model, the XFEM bypasses the need for introducing fictitious materials and allows directly describing the constitutive behavior of each phase by appropriate models. In this study, we build upon our recent work on combining the LSM and the XFEM [13–15] and expand this approach onto the optimization of PACs.

To accurately model the response of PACs, in general, large deformations and a complex thermomechanical constitutive behavior need to be considered, leading to a nonlinear nonconservative multiphysics model. Constitutive models for accurately describing the SM cycle of amorphous SMPs are presented, for example, in Refs. [1], [5], [16], and [17]. Considering these models within the topology optimization process significantly increases the algorithmic complexity and the computational cost. To create an analysis approach more amenable for design in terms of computational speed and ease of implementation but with acceptable fidelity, we introduce two simplified mechanical models. Both models assume small deformations but differ in approximating the residual strains due to the programing loads. We study the accuracy of these PAC models against experimental results and compare the optimization results obtained using these models.

The primary goal of this study is introducing a LSM–XFEM topology optimization approach for the design of PACs. The main components of this approach are illustrated in Fig. 1. The geometry of the active and passive material phases is described by a discretized LSF; its parameters are altered in the design process via a gradient-based optimization scheme. The response of a given material layout and the design sensitivities are predicted by an XFEM discretization of a PAC model. The XFEM analysis uses the same mesh used for discretizing the LSF. The triangular elements shown in Fig. 1 are only used for integrating the weak form of the governing equations over the individual phases and do not require to remesh the PAC design as the geometry evolves in the optimization process. Upon convergence of the design, the XFEM mesh used for integration is extracted along the material interfaces and the outer boundaries, saved in an appropriate file format, and loaded into the printer interface. To distinguish between the individual phases, a separate surface mesh needs to be extracted for the active and passive phases.

The underlying theories, algorithmic details, and numerical examples of the proposed optimization method are described in the remainder of this paper that is organized as follows: The physical phenomena dominating the response of PACs are outlined and the simplified mechanical models are described in Sec. 2. The LSM–XFEM design optimization approach is summarized in Sec. 3. The PAC models and the LSM–XFEM framework are validated against experimental results in Sec. 4. The optimization approach is studied via numerical examples in Sec. 5. The main findings are discussed in Sec. 6.

To provide insight into the nature of the design problem, we describe the fabrication of PACs and summarize the physical phenomena that need to be considered in the PAC design. The formulations of simplified multiphysics models approximating the PAC response upon activation are presented.

The PACs considered here are fabricated by 3D printing a heterogeneous layout of two materials, a glassy SMP and an elastomer [1,18]. The activated shape of the PAC is programed by a sequence of thermomechanical processing steps depicted in Fig. 2. After the PAC is printed, it is heated to a temperature T_H that is above the glass transition temperature T_g of the SMP phase, i.e., T_H > T_g. The PAC is then deformed by applying prescribed displacements. We refer to this configuration as programing stage (i). While the loading is maintained, the PAC is then cooled to a temperature T_L that is below T_g, i.e., T_L < T_g. Finally, the loading is removed and the PAC deforms. This configuration is referred to as programed stage (ii). Upon heating the PAC above T_g it recovers its initial shape. Note that the topology of the phases in the composite is arbitrary; Fig. 2 shows just one possible instantiation: SMP fibers in an elastomeric matrix. While the PAC may be initially planar, as shown in Fig. 2, the 3D printing process as well as our computational design methodology allows the fabrication of PACs with arbitrary 3D geometries.

The deformation of the PAC upon releasing the prescribed loads is due to residual stresses generated during the thermomechanical processing steps. At T_H, both the SMP and elastomeric phase are in a soft, rubbery state with similar elastic properties. As the temperature is lowered and the SMP phase undergoes glass transition, the stresses in the SMP phase are partially relaxed and, at the same time, the material stiffness of the SMP phase increases [17]. The level of stress relaxation in the SMP phase is characterized by the fixity, f, with f = 1 representing a complete stress relaxation. The phase change in the SMP phase causes residual stresses in both the SMP and matrix phases as the external loads are removed. In the example of Fig. 2, the SMP phase is in compression and the matrix phase in tension upon unloading. The shape of the PAC in stage (ii) depends on the initial shape, the layout of the individual material phases, and the programing loads. Here, we limit the complexity of the design problem by considering only initially planar configurations and simple loading scenarios that are not subject to optimization.

The particular fabrication technology considered in the design studies in Sec. 5 is the multimaterial polymer printing process of an Objet Connex 260 printer (Stratasys, Edina, MN). The PACs are printed by depositing droplets of polymeric ink at 70 °C, wiping them into polymeric film with a roller, and then UV photopolymerizing the film by multiple passes of UV lamps. Simultaneous deposition of two distinct polymers in preset combinations enables creation of parts with graded material properties. In a recent study [3], we exploited this design freedom to create PACs, but to simplify matters here we only consider two distinct material compositions: Tango black + for the passive phase and Shore95 for the active phase. The material properties of these materials are given in Table 1. While both materials exhibit a SM behavior, operating the PAC above −5 °C ensures that the passive phase remains in a soft rubbery state. Here, we consider material properties for a SM cycle between T_L = 15 °C and T_H$≥$ 60 °C.

To predict the shape of the PAC in the programed stage (ii), we study two simplified PAC models. Both models account for stress relaxation and the resulting residual stresses via eigenstrains that can vary with position. These models lump the complex nonlinear thermomechanical phenomena occurring during a thermomechanical cycle into an effective inelastic strain, i.e., eigenstrain. Both PAC models assume infinitesimal strains and a linear material behavior, but they differ in the approximation of the eigenstrains in the active phase. The more accurate model predicts the eigenstrains by resolving the strain field in the PAC at stage (i) due to the programing loads at T_H. The simpler model assumes a uniform eigenstrain in the active phase which is predicted analytically. Subsequently, we first describe the more accurate model which is then simplified to the model with uniform eigenstrains. Both models are compared against experimental data in Sec. 4.

We consider the response of the PAC at stage (i) subject to programing loads and at stage (ii) after the programing loads are removed. The effects of the phase transition of the active phase are modeled via eigenstrains, which are computed in stage (i) and applied in stage (ii), and a change in Young's modulus. The residuals of the static equilibrium equations in both stages are written in weak form as follows:

with

where $(·¯)$ denotes quantities in stage (i) and $(·∧)$ in stage (ii). Passive and active phases are denoted with superscripts “1” and “2” as shown in Fig. 3. The displacement vector, the strain tensor, and the stress tensor are denoted by $uip,ɛijp$⁠, and $σijp$⁠, respectively. External forces, $Fip$⁠, act along the boundary $ΓTp$⁠, and prescribed displacements, $Uip$⁠, are applied at the boundary $ΓUp$⁠. Note the boundary conditions typically differ between stages (i) and (ii). The material behavior at the two stages is described by the following constitutive equations:

where $Cijklp$⁠, f^p, and α^p denote the elastic tensor, the strain fixity, and the coefficient of thermal expansion, respectively. The Young's modulus of the active phase increases as the material transitions from a rubbery state in stage (i) to a glassy state in stage (ii). The elastic tensor of the passive phase is assumed constant. The strains of the programing state (i) are introduced as eigenstrains in Eq. (5). For the passive phase, the fixity is assumed negligible, i.e., $f1=0.0$⁠. In the following, we ignore thermal strains as the mismatch in α^p is insignificant for materials considered in this study; see Table 1.

Subsequently, we refer to the model described above as “two-stage” model. This model can be further simplified by approximating the strains in the programing stage (i) by a uniform strain distribution. Assuming a homogeneous material and for simple geometries and programing loads, $ɛ¯ijp$ can be approximated analytically via kinematic compatibility considerations. For example, in the configuration depicted in Fig. 3, the strain along the bar axis is simplified to $ɛ¯11p=U¯1/L$⁠, where L is the length of the bar. Assuming nearly incompressible material, the transverse eigenstrains are taken as $ɛ¯22p=ɛ¯33p=-0.5ɛ¯11p$⁠. All other eigenstrains are assumed to be zero. As before, the fixity in the passive phase is assumed to vanish. Note that this simplification requires only modeling the static equilibrium of stage (ii) and solving Eq. (2). We refer to this model as “single-stage” model.

To optimize the spatial arrangement of passive and active phases, we adopt a topology optimization approach that describes the material layout by a LSM and predicts the PAC response by an XFEM discretization of the mechanical models described above. To this end, we adopt our LSM–XFEM approach previously studied for problems in fluid and structural mechanics [13–15,19] and study it for the design of PACs. The main concepts of the LSM–XFEM approach are outlined in this section. For a detailed discussion of the relationships between our LSM–XFEM approach and other structural optimization methods using level sets and immersed boundary techniques, the reader is referred to Ref. [14].

The geometry of phases 1 and 2 is described with a LSF $φ(xi)$⁠, where x_i denotes the spatial coordinate. Negative values, $φ(xi)<0$⁠, define the domain $Ω1$ occupied by phase 1, positive values, $φ(xi)>0$⁠, the domain $Ω2$ occupied by phase 2, and the zero level set contour, $φ(xi)=0$⁠, defines the phase boundary $Γ1,2$⁠. The LSF can be parametrized such that it describes geometric primitives and the optimization variables define location, orientation, and dimensions of the primitives. This approach is used to describe the layout of active fibers in the validation study of Sec. 4. Alternatively, to increase the design freedom and to allow for the evolution of complex geometries in the optimization process, the LSF is discretized by a finite-element mesh and the nodal level set values are defined as functions of the optimization variables. The resulting parameter optimization problem is solved with a nonlinear programing method. The latter approach is used in the design optimization studies of Sec. 5.

To accelerate convergence of the optimization process, a linear filter is applied to the nodal level set values [13]. Assigning one optimization variable s_j to each node of the finite-element mesh, the nodal LSF values, $φi$⁠, are defined as follows:

where N_n denotes the number of nodes in the finite-element mesh, r_s is the smoothing radius, and $|xi-xj|$ is the distance between nodes i and j. The level set filter (6) widens the zone of influence of the optimization variables on the level set field and thus enhances the convergence of the optimization process [13].

The main advantage of the LSM is the crisp definition of the phase boundaries, $Γ1,2$⁠. However, the LSM is in essence driven by shape sensitivities which do not promote the creation of new inclusions of one phase within volumes occupied by the other phase. This along with the nonconvexity of topology optimization problems often leads to a strong dependence of the optimization results on the initial design. To mitigate this issue, the optimization process can be augmented by seeding new inclusions, for example, via topological derivatives [20–22]. In this study, we adopt a different strategy and address the issue by starting the optimization process from a dense array of inclusions.

The structural response of the PAC is predicted by the XFEM. This discretization method retains the crisp description of the geometry defined by the LSM and operates on fixed meshes which do not change in the optimization process. The latter enhances the robustness and efficiency of the structural analysis.

The basic concept of the XFEM is illustrated in Fig. 4. We consider a finite-element mesh that is not aligned with the interfaces of a given material distribution described by a LSF. The two elements connected to node 3 are intersected by the zero level set isocontour, i.e., the material interfaces. To capture the discontinuities in the strain and stress fields across the phase boundaries, the standard finite-element approximation spaces are enriched with additional shape functions [12]. Here, we use a Heaviside enrichment approach [23] and approximate the response in each phase separately by standard finite-element shape functions. Solid lines mark the domains where the displacements are interpolated by the shape function associated with a particular phase and dashed lines indicate that the displacements are not interpolated by this shape function. As this approximation allows for discontinuities across the phase boundaries, the continuity of the displacement field needs to be enforced explicitly.

To consistently approximate the displacement field in the presence of small geometric features and to remedy artificial coupling between disconnected inclusions, we adopted a generalized formulation of the Heaviside enrichment strategy [14]. The displacement fields of the single-stage and two-stage PAC models defined by Eqs. (1) and (2) are approximated as follows:

with H being the Heaviside function

where N_z is the local shape function, $Nnodese$ is the number of nodes of the eth element, and $ui,mq,z$ is the displacement degree-of-freedom (DOF) at node z for phase $q=[1,2]$ in the ith direction. The enrichment level is m, with M being the maximum number of enrichment levels. The Heaviside function, H, turns on/off two sets of shape functions associated with the phases 1 and 2. The Kronecker delta, $δrmz$⁠, selects the enrichment level, r, such that the solution at a point, x, is interpolated only by a single-shape function defined at node z. For each phase, multiple enrichment levels, i.e., sets of shape functions, are necessary if the DOFs interpolate the solution in multiple, physically disconnected regions; see Refs. [14], [24], and [25]. This generalization prevents spurious coupling and load transfer between disconnected regions of the same phase [14].

To enforce continuity of the displacements across phase boundaries, we augment the weak forms of the governing equations (1) and (2) by a stabilized version of the Lagrange multiplier method

with

and

where $λ$ is the Lagrange multiplier, and $δλ$ is the associated test function. Increasing the weighting factor, γ, improves the accuracy of enforcing the continuity of the displacement field across the phase boundary, at the cost of higher condition numbers and reduced numerical stability. The XFEM may also suffer from ill-conditioning as the level set field creates small intersections within an intersected element. We cure this issue by a geometric preconditioning scheme [26].

While the formulation of the proposed LSM–XFEM approach admits a broad range of objectives and multiple design constraints; here, we focus on one particular class of problems: designing the microstructure of a PAC such that the composite deforms into a target shape at stage (ii). To this end, we minimize the square of the difference between the displacements in the programed stage (ii), $u∧i$⁠, and a target deformation $uitarget$ integrated over the surface $Γtarget$ where the deformations are monitored. This yields the following contribution to the objective function:

where w_i denotes the weighting factors for the components of the displacement vector. To prevent the formation of small geometric features and irregular shapes, we augment the objective (13) with the following regularization term that penalizes the perimeter of the phase boundary:

where γ_per denotes the perimeter penalty factor. The value of the perimeter penalty factor is chosen to smooth the surface of the phase boundary without significantly affecting the shape matching objective $Ftarget$⁠. The choice of γ_per is problem dependent and typically needs to be determined iteratively. Note that while penalizing the perimeter often leads to satisfactory results for a broad range of design problems, it does not allow for explicitly controlling the minimum feature size [11,14].

The formulation of the optimization problem can be summarized as follows:

where N_s is the number of optimization variables, and s_min and s_max denote the lower and upper bounds on the optimization variables, respectively. The box-constrained optimization problem (15) is studied in Sec. 5.

The PAC models introduced above are built upon several simplifying assumptions. To investigate the ability of these models to capture the fundamental characteristics of the PAC response, we consider the simple configuration of Fig. 5, which was studied experimentally by Ge et al. [1]. The active fibers are made of Shore95 and exhibit a SM behavior; the passive phase is Tango black+ . In the programing stage (i), the composite plate is subject to prescribed displacements $±0.5U¯1$ at the front and rear faces. In the programed stage (ii), it is clamped at its center to suppress rigid body motion. In the case of the single-stage model, the eigenstrain is assumed to be $ɛ¯ij=diag[U¯1/L,-0.5U¯1/L,-0.5U¯1/L]$⁠, which is consistent with the programing displacements and incompressible material. The length of the plate in x₁-direction is denoted by L.

The arrangement of uniformly spaced fibers is described by the following LSF:

with r_f, $ϕ$⁠, and Δ_f being the fiber radius, the fiber orientation, and the distance between fibers. The relative location of the ith fiber is determined by the parameter $ci=[-6:1:6]$⁠. The value $di(x1,x2)$ defines the distance to the centerline of the ith fiber in the plane $x3=0$⁠. Finding the minimum distance $di2(x1,x2)+x32$ allows the definition of an array of uniformly spaced fibers by a single equation. Figure 5 shows the discretized LSF and the XFEM model of active and passive phases for a fiber orientation of $ϕ=30deg$⁠. The nodal level set values are computed by Eq. (16) and interpolated by trilinear shape functions within the finite elements.

The parameters of the XFEM model and the material properties are given in Tables 2 and 3, respectively. The computational mesh consists of 33,614 elements and 39,600 nodes. For the fiber orientation $ϕ=30deg$⁠, 10,976 elements are intersected resulting in a XFEM model with 370,154 DOFs.

The response of the PAC is studied for seven different fiber orientations. The experimental results of Ge et al. [1] and the XFEM predictions are compared in Fig. 6. The fixity of the composite plate is measured by the elongation of the plate in the programed stage (ii) relative to prescribed elongation $U¯1$ at the programing stage (i). The fixity parameter of the active phase is calibrated by considering the experimental results for a plate made of the active phase only; see the dotted line in Fig. 6. Comparing the experimentally measured and numerical results shows that the two-stage model captures well the PAC response across a wide range of fiber orientations. However, the single-stage model leads only to acceptable results for fiber angles up to $ϕ=20deg$⁠, because the assumption for the eigenstrain distribution is valid only for small differences between training load and fiber directions.

Given the complexity of the material behavior and the uncertainty of the experimental results [1], the two-stage model and, to a lesser extent, the single-stage model show satisfactory agreement with the experimental measurements for the example considered above. This motivates the application of these PAC models to the design optimization problems presented in Sec. 5. However, additional validation problems need to be studied to increase the confidence in the proposed PAC models. Future studies should also include problems that exhibit large out-of-plane deformations and more complex inclusion geometries.

We study the proposed LSM–XFEM approach with four design problems, using the single-stage and two-stage PAC models. In the first three examples, the target shape promotes a plate-bending type response of increasing complexity. The fourth example explores the ability to generate localized deformations. Measurements for 3D-printed realizations of the optimized designs are compared against numerical predictions of the programed shape at stage (ii). The design studies illustrate the ability of the proposed optimization approach to find nontrivial designs. The comparison against experimental measurements motivates future research on improving the accuracy of the PAC model.

For all examples, we study the formulation of the optimization problem (15) considering different target shapes and PAC models. The geometry of the samples and boundary conditions in the programing stage (i) are depicted in Figs. 7 and 8 for the plate bending and for the localized deformation problems, respectively. The dimensions and mesh sizes are listed in Table 4. The difference between the target and programed shape (13) is measured at the upper surface for the first three examples and at the upper and lower surfaces for the fourth problem. The analytical expressions for the target displacements are given in Table 5 and are defined with respect to a coordinate system located at the center of the design domain.

In the programing stage (i), the samples are stretched by 6.0% by applying prescribed displacements on front and rear faces. In the programed stage (ii), the samples are supported at their center to prevent rigid body motion. The material parameters are given in Table 3. In all examples, the design domain is discretized by eight-node, trilinear finite elements. The parameters of the XFEM meshes are given in Table 4. The number of DOFs in the XFEM model depends on the intersection configurations and varies between 179,000 and 300,000 DOFs for the first three examples and 350,000 and 422,000 for the fourth example. The linear systems of equations resulting from the discretized PAC models are solved by a Generalized Minimal Residual Method (GMRES) solver, preconditioned with an incomplete LU decomposition [27].

The LSF is discretized by the same mesh used for the PAC analysis. The nodal level set values are defined as functions of the optimization variables, using the linear filter (6). The smoothing radius is set to $rs=1.6 mm$ for the first three examples and to $rs=0.64 mm$ for the fourth example. The resulting parameter optimization problems are solved by the Globally Convergent Method of Moving Asymptotes (GCMMA) of Ref. [28] without subcycling. The parameters for the optimization problem (15) and for the GCMMA are given in Table 6. The design sensitivities of the objective and constraints are determined by the adjoint method. The partial derivatives of the state equations and objective function with respect to the state variables are evaluated using analytically differentiated formulations. The partial derivatives of the objective, constraints, and element residuals with respect to the optimization variables are calculated by a finite-difference scheme. Note that the computational cost of the finite-difference scheme is insignificant, as the partial derivatives vanish for all but the optimization variables affecting the shape of the phase boundaries. The adjoint problems are also solved by a gmres solver, preconditioned with an incomplete LU decomposition.

The designs are initialized with uniform arrays of cuboids of either active or passive material, surrounded by material of the other phase. The arrays are defined by the following LSF:

where a, b, and c denote the dimension the cuboids in x₁-, x₂-, and x₃-directions, respectively. The coordinates of the center of the individual cuboids are $[x1,i,x2,i,x3,i]$⁠. The roundness of the cuboids is determined by the exponent n which is set to n = 100.0. In the first three examples, the optimization processes are started from 2D arrays of cuboids placed at the middle plane of the plate. In the fourth example, a 3D array of cuboids is used. The material phase and the dimensions of the cuboids, as well as the size of the array of initial inclusions, are given with each example.

The optimized designs are 3D printed by exporting the triangulated phase boundaries of the XFEM models for each phase separately, using the stereolithography (STL) file format. To this end, we extract the XFEM mesh used for integrating the weak form of the governing equations. Note that this triangulation is aligned with the phase boundaries and thus each STL file defines the volume of one distinct phase; see Sec. 3.2. The STL files are loaded into the software interface of an Objet Connex 260 printer and a material type is assigned to each volume associated with a STL file. The printer software optimizes the material layout along the phase boundaries for improved bonding and material strength. The printed PAC sample is subject to the thermomechanical programing cycle described above and the programed shape in stage (ii) is measured. Note that the experimental results reported here mainly serve as a proof of concept but lack the accuracy needed for rigorous quantitative comparisons between predictions and experimental measurements.

In this example, we seek to find the geometry of active inclusions such that the programed PAC plate assumes a parabolic shape in stage (ii); see Table 5. The target shape corresponds to a maximum displacement magnitude of $|u3target|=9.6 mm$ along the edges $x1=±40 mm$ ; the target displacement along the symmetry line $x1=0.0 mm$ is zero. We start the optimization process with a 6 × 4 array of cuboids of active material surrounded by passive material; see Eq. (18). The dimensions of the cuboids are $a=4.0 mm,b=2.4 mm$⁠, and $c=0.8 mm$⁠. We enforce a symmetric design about the x₁–x₃ and x₂–x₃ planes via the definition of the nodal level set values as functions of the optimization variables. The XFEM analysis is performed on the full mesh. We optimize the material layout with the single-stage and two-stage models and study the influence of the PAC models on the optimized designs. Measurements for a 3D-printed realization of the optimized two-stage design are compared against numerical predictions of the programed shape in stage (ii).

The optimized material layouts are depicted in Fig. 9. The deformations of the optimized designs at the programed stage (ii) are shown in Fig. 10. The shape mismatches of the initial and optimized designs are reported in Table 7. The evolutions of the shape mismatch (13) and the perimeter in the optimization process are given in Fig. 11.

For both PAC models, the designs converge in less than 250 iterations and result in a substantial reduction of the shape mismatch (13). For both models, the optimization process leads to a fiberlike arrangement of the active phase. The geometry of the active inclusions is rather complex but can be realized by 3D printing without manufacturing complications. The design optimized with the two-stage model leads to wider, curved fibers that are connected via bulky sections at the ends of the plate. The latter feature stiffens the plate and reduces cylindrical bending about the x₁-axis. Therefore, the design optimized with the two-stage model utilizes a larger amount of active material.

The evolutions of the shape mismatch and the perimeter in the optimization process are smooth. Toward the end of the optimization process mainly the perimeter is reduced while the shape mismatch error remains nearly constant. A similar behavior can be observed for the other examples discussed subsequently. In average, one optimization step with the two-stage model, including a forward analysis and an adjoint sensitivity analysis, requires 12.1 min using four cores of a state-of-the-desktop computer. For the single-stage model, the time per iteration step drops to 2.9 min as this model leads to less DOFs in the XFEM analysis and computing the derivatives of the elemental residual is less costly. Note that the in-house code used for this study is not optimized for speed. The computational cost for the remaining examples is comparable to the one reported above.

To further illustrate the influence of the PAC models, we analyze the design optimized with the single-stage PAC model by the two-stage model and vice versa. The shape mismatch (13) of the two-stage optimized design is $826.9 mm4$ when analyzed with the single-stage PAC model; and the one of the single-stage optimized design is $819.0 mm4$ when analyzed with the two-stage PAC model. The mismatch of the single-stage optimized design is more than seven times that of the two-stage optimized design, when both designs are analyzed with the more accurate two-stage model; see Table 7. This comparison highlights the need for accounting for the inhomogeneous strain fields in the programing stage (i) when optimizing the PAC design.

Finally, we compare the XFEM predictions for the design optimized with the two-stage PAC model against experimental measurements of the printed design which is shown in Fig. 12. Note for visualization purposes, we print the sample with a transparent matrix material and a dark SMP material. However, the mechanical test is performed with the materials specified above. Qualitatively, the printed design assumes a parabolic shape without significant cylindrical bending about the x₁-axis. The maximum measured deflection is about 18% larger than the XFEM prediction. This discrepancy can be attributed to the simplifications of the two-stage PAC model and to errors of our desktop experiments. Considering the magnitude of the displacements relative to the plate thickness, more accurate PAC models that account for finite strains and large displacements should be used in the optimization process.

We repeat the design study described in Sec. 5.1 but consider a cosine wave as target shape; see Table 5. Here, the programed shape includes concave and convex sections which are expected to require a more complex layout of active material. To accommodate the increase of complexity in the target shape, the magnitude of the maximum target displacement is reduced to $4.8 mm$⁠. Again, we enforce a symmetric design about the x₁–x₃ and x₂–x₃ planes via the definition of the nodal level set values as functions of the optimization variables. The XFEM analysis is performed on the full mesh. We start the optimization process with a 6 × 2 array of cuboids of active material surrounded by passive material. The dimensions of the cuboids are $a=4.0 mm,b=2.4 mm$⁠, and $c=0.8 mm$⁠. As in the previous example, we optimize the material layout with the single-stage and two-stage PAC models.

The optimized material layouts are depicted in Fig. 13. The deformations of the optimized designs at the programed stage (ii) are shown in Fig. 14. The shape mismatches of the initial and optimized designs are reported in Table 8. As in the previous example, the designs converge quickly and result in a substantial reduction of the shape mismatch (13). Both PAC models lead to fiberlike arrangements of the active phase where the fibers are closely aligned with the direction of the programing loads. The two-stage optimized design features sophisticated arrangements of active material near the ends of the plates, which facilitate the transition between the convex and concave bending and counteract cylindrical bending about the x₁-axis. As in the previous example, the design optimized with the single-stage model utilizes a smaller amount of active material.

The impact of the differences in the geometry between the optimized designs is illustrated by analyzing the design optimized with the single-stage PAC model by the two-stage model and vice versa. The shape mismatch (13) of the two-stage optimized design is $656.8 mm4$ and for the single-stage optimized design $606.0 mm4$⁠. As in the previous example, the increase in shape mismatch is significant. The shape mismatch of the single-stage design is more than eight times that of the two-stage optimized design, when both designs are analyzed with the more accurate two-stage model.

In Fig. 15, we show the 3D-printed design optimized with the two-stage PAC model. As previously, we print the sample with a transparent matrix material and a dark SMP material only for visualization purposes, but we perform the mechanical test with the materials specified in Table 3. Qualitatively, the printed design assumes the target shape. However, the transition from the concave center section to the concave outer sections is less pronounced than predicted by the XFEM analysis. The maximum measured deflection is about 75% larger than the XFEM prediction, further emphasizing the need to consider more accurate PAC models in the optimization process as well as improving the accuracy of the experiments.

To further increase the complexity of the target shape, we consider on overlay of a parabolic bending about the x₂ and a twist about the x₁-axis as target shape; see Table 5. The magnitude of the maximum deflection in x₃ direction, which occurs at the corner points, is chosen to be $4.8 mm$⁠. Here, we start the optimization process with a 6 × 4 array of cuboids of passive material embedded in active material. The dimensions of the cuboids are $a=4.0 mm,b=2.4 mm$⁠, and $c=0.8 mm$⁠. We enforce a symmetric design about the x₂–x₃ plane. The XFEM analysis is performed on the full mesh. As before, we optimize the material layout with the single-stage and two-stage models.

The optimized material layouts are depicted in Fig. 16. The deformations of the optimized designs at the programed stage (ii) are shown in Fig. 17. The shape mismatches of the initial and optimized designs are reported in Table 9. In contrast to the previous examples, the designs converge slower but eventually result in a substantial reduction of the shape mismatch (13). Furthermore, the active phase is not arranged in long fibers that are aligned with the direction of the programing loads. Instead, the active phase forms fibers oriented about $±45deg$ with respect to the x₁-axis. The twist deformation is facilitated by an offset of about 90 deg between the top and bottom fibers. While both PAC models lead to similar arrangements of the active phase, the geometry of the active phase differs in particular along the outer edges and the plate center.

We analyze the performance of the two-stage optimized design with the single-stage model and vice versa. The shape mismatch (13) of the two-stage optimized design increases to $725.9 mm4$ and the one of the single-stage optimized design to $589.83 mm4$⁠. The shape mismatch of the single-stage optimized design is more than three times that of the shape matching error of the two-stage optimized design, when both designs are analyzed by the two-stage model. While this difference is less dramatic than in the previous examples, it still motivates the use of PAC models that resolve the strains at the programing stage.

A 3D-printed sample of the design optimized with the two-stage PAC model is shown in Fig. 18. As previously, we print the sample with a transparent matrix material and a dark SMP material only for visualization purposes, but we perform the mechanical test with the materials specified in Table 3. Qualitatively, the printed design assumes the target shape. The measured deflections of the corner points are about 30% larger than predicted by the XFEM.

In the final example, we study the ability of the proposed optimization approach to find a PAC design that yields a localized deformation, defined by a single parabolic bulge; see Table 5. The amplitude of the desired bulge is $0.1 mm$⁠. Note for this problem, the target shape is defined on both the upper and lower surface of the design space; see Fig. 8. The design space is restricted to a subdomain of the composite plate, which is shorter and thicker than the one considered in the previous examples. To reduce the computational cost, we exploit the symmetry of the design problem and consider only one-eighth of the design domain in the optimization process. We initialize the LSF in the one-eighth model with a $4×4×4$ array of cuboids of active phase. The dimensions of the cuboids are $a=0.4 mm,b=0.4 mm$⁠, and $c=0.2 mm$⁠. Again, we optimize the material layout with the single-stage and two-stage models.

The optimized material layouts are shown in Fig. 19. The deformations of the optimized designs at the programed stage (ii) are given in Fig. 20. The shape mismatches of the initial and optimized designs are reported in Table 10. These results suggest that designing PACs to achieve a localized deformation is more challenging. The target shape is not matched well by neither designs. The maximum amplitudes of the bulge for the designs optimized with the single-stage and two-stage models reach only 68.8% and 38.3% of the target value. Nevertheless, the sophisticated spatial arrangements of the active phase demonstrate the ability of the proposed optimization method to resolve geometric details. For both PAC models, the active phase is arranged to create a compliant mechanism that pulls the center of the top and bottom surfaces inwards. However, the stiffness of the passive phase prevents large deformations. The low volume fractions of the optimized designs suggest that larger deformations cannot be obtained by just increasing the amount of active phase.

Analyzing the two-stage optimized design with the single-stage PAC model and vice versa yields a shape mismatch of $0.17 mm4$ for the two-stage design and $0.36 mm4$ for the single-stage design. The shape matching error of the single-stage optimized design is 44% larger than one of the two-stage optimized design, when predicting the shape mismatch with the two-stage PAC model. This result suggests that even when small localized deformations are desired the PAC model has a significant influence on the optimized design and an accurate PAC model should be used for design optimization.

In Fig. 21, we show the 3D-printed design optimized with the two-stage PAC model. Here, we print the sample with a transparent matrix material and a dark SMP material. As our experimental instrumentation lacks the ability to sufficiently resolve submillimeter deformations, we do not report on experimental measurements for this problem.

This paper has introduced a topology optimization approach for the design of 3D PAC s. A LSM was used to describe the geometry of boundaries between an active SMP phase and a passive elastomeric matrix phase. To avoid the complexity of previously proposed PAC models, two simplified models were introduced to predict the PAC response, both assuming small strains and a linear material behavior. The single-stage model approximates the residual strains due to the programing loads by a uniform eigenstrain distribution. The two-stage model predicts the strains in the programing stage and utilizes them as residual strains to describe the deformed shape of the programed PAC. A generalized formulation of the XFEM was adopted to discretize the PAC models. Comparing XFEM predictions for both PAC models against experimental results for a PAC plate with active fibers suggested that the two-stage model is valid for a wide range of fiber orientations. However, the single-stage model is only applicable to configurations where the active inclusions are closely aligned with the direction of the programing loads.

The design optimization studies focused on finding the microstructural layout such that a PAC plate assumes a target shape. The LSF was discretized by finite elements and the nodal level set values were defined as explicit functions of the optimization variables using a linear filter. The objective was composed of a measure of the shape mismatch and a penalty on the perimeter that regularizes the problem. The weighting of the penalty term was found problem dependent and determined iteratively.

The optimization examples included three problems where the target shape corresponds to a plate-bending/twisting type deformation and one problem where a localized deformation was desired. The numerical results demonstrated that the proposed optimization approach provides a crisp and well-resolved description of the phase boundaries. Surface meshes of the phase boundaries can be easily extracted from the XFEM model and directly exported into the software interface of 3D printers without the need for any form of postprocessing.

In the first three problems, the target shapes could be matched with acceptable accuracy. Both PAC models resulted in conceptually similar designs with sophisticated spatial arrangements of the active phase. The results for the fourth example suggest that creating localized deformations is a challenging design problem even if the amplitude of the target design is small. Nevertheless, these examples demonstrated the ability of the proposed optimization method to find nontrivial and nonintuitive designs. Analyzing the designs optimized with the single-stage PAC model by the more accurate two-stage PAC model showed that the shape mismatch is significantly larger than the one of the designs optimized with the two-stage model. This observation suggests that the accuracy of the PAC model has a noticeable influence on the optimization results and therefore an accurate PAC model should be used for topology optimization.

The comparison of the XFEM predictions against experimental measurements for the optimized designs of the first three examples suggested that the proposed two-stage PAC model may be appropriate for design purposes but insufficient to accurately predict the PAC response in the presence of large deformations. In future studies, more accurate but still computationally tractable PAC models will be researched and integrated into the proposed optimization framework. Furthermore, the perimeter constraint was found to be insufficient to explicitly control the local feature size of the active phase. Thus, methods will be explored to constrain the minimum feature size and to account for the resolution limits of 3D printers.

The first four authors acknowledge the support of the Air Force Office of Scientific Research under Contract No. FA9550-13-1-0088 and the National Science Foundation under Grant No. EFRI-ODISSEI 1240374. Maute and Qi further acknowledge the support of the National Science Foundation under the collaborative research Grant No. CMMI 1463287/1462894. Ding and Dunn acknowledge the support of the SUTD Digital Manufacturing and Design (DManD) Centre supported by the National Research Foundation of Singapore. The opinions and conclusions presented in this paper are those of the authors and do not necessarily reflect the views of the sponsoring organization.