Programming temporal morphing of self-actuated shells

Advances in shape-morphing materials, such as hydrogels, shape-memory polymers and light-responsive polymers have enabled prescribing self-directed deformations of initially flat geometries. However, most proposed solutions evolve towards a target geometry without considering time-dependent actuation paths. To achieve more complex geometries and avoid self-collisions, it is critical to encode a spatial and temporal shape evolution within the initially flat shell. Recent realizations of time-dependent morphing are limited to the actuation of few, discrete hinges and cannot form doubly curved surfaces. Here, we demonstrate a method for encoding temporal shape evolution in architected shells that assume complex shapes and doubly curved geometries. The shells are non-periodic tessellations of pre-stressed contractile unit cells that soften in water at rates prescribed locally by mesostructure geometry. The ensuing midplane contraction is coupled to the formation of encoded curvatures. We propose an inverse design tool based on a data-driven model for unit cells’ temporal responses.


Supplementary Note 1: Shell design
Our shells have three layers. Two 3D-printed non-uniform tessellations form the outer layers. These are glued to either side of a pre-stretched latex membrane. This section describes the design of the flat-printed geometries.

Single unit cell design
The shells' outer layers are tessellations of the structures shown in the Supplementary Figure 1. Each of these unit cells is bounded by four cylinders of diameter 4 mm and height 2.3 mm, called bases. The bases are connected to their neighbors by two symmetric V-shaped spring elements (called brackets) that form an angle of 37 • relative to the central axis and have a height of 1 mm. At room temperature, brackets are sufficiently stiff to prevent finite deformations due to compression by the pre-stretched elastic membrane that constitutes the shells' mid-planes. The brackets soften when placed in hot water, inducing in-plane contraction. This contraction occurs until the bumpers attached to the bases in the space between the brackets collide. This collision occurs in each unit cell once the target in-plane deformation for that unit is reached. Local curvatures are programmed in a unit cell by setting different bumper lengths for the opposite outer layers. To facilitate aligning the two outer layers with respect to each other during shell fabrication, cylindrical holes are subtracted from several bases (see the Supplementary Note 5 for more information). The feature dimensions were chosen for the following reasons. Bases interface the membrane to the brackets, so they must have a sufficiently large gluing area to be reliably connected to the membrane, but should be small enough to allow large curvatures in decimeter-scaled specimens. The bracket shape is designed to reduce both in-plane and out-of-plane shearing. In-plane shearing is prevented due to the large (37 • ) angle between the bracket and the central axis. Out-of-plane shearing is prevented due to the rectangular shape of the bracket section. The 1 mm bracket height is always larger than its thickness (h ≤ 0.65 mm), which makes in-plane bending energetically favorable.

Tessellating unit cells
The design pipeline starts with a user-provided target surface (Supplementary Figure 2), which is isotropically remeshed into a triangular mesh T with a target number of vertices N [1]. Each pair of adjacent triangles of this mesh represents one unit cell with two bases placed on opposite faces of each triangle, centered on its barycenter. Thus, the whole mesh T serves as a stencil for our structure's final state, when deformation has completed. It is scaled to ensure there is a minimal bumper length (0.3 mm) to limit bracket deformations in the final state. The actual size of the shell depends on the number of vertices in T . In order to generate a shell of dimensions close to the input surface, we aim to minimize the required scaling of the stencil. First, we find the best fitting number of vertices of the stencil. It can be coarsely estimated as N ≈ 0.38A/l 2 avg where A is total surface area and l avg = 7 mm is the average unit cell length (using heuristic knowledge that the number of triangles is twice the number of vertices, assuming triangles are close to regular, and approximating unit cell lengths as twice the radius of circumscribed circles). Then, N can be varied to find a value leading to minimal stencil scaling.
Given the layout of bases on the stencil, the bumpers are first constructed as boxes that are aligned along the lines connecting the centers of the triangles and projected onto the corresponding stencil triangle. Matching pairs of boxes are trimmed by the bisector plane between the triangles, defining the interface between neighboring bumpers.
The flat arrangement of bases is then constructed using a minimal distortion conformal map [2] from T to a resulting 2D mesh, F. We exploit conformal flattening since it circumvents shearing which would result in undesired shear forces in the membrane as it contracts. The mesh F serves as a stencil for the structure's initial, printed state. We relocate the bases by translation from triangle centers in the final state's stencil to triangle centers in the initial state's stencil. The bases are rotated to align one of its bumper axes in the direction of the corresponding neighbor's center since it is generally not possible to perfectly align all of them. Then, F is scaled so that the relocated bases with bumpers do not intersect with their neighbors, moreover specifying the minimal gap between them necessary for fabrication (0.1 mm).
Once the flat layout of bases with bumpers is complete, V-shaped brackets that bridge the gaps between bases are generated. This produces a flat structure that fully encodes the target geometry. Conformal flattening may produce overlaps of unit cells, or some of the unit cells may be longer than our upper bound (9.5 mm), making the target shape impossible to replicate. Large unit cell lengths may occur if the target surface has regions with high curvatures or that require substantial stretching during conformal flattening. It is in certain cases possible to resolve the latter issue by placing cuts across the surface [3] or by editing the input geometry. Morphing times are programmed by configuring the thicknesses of the brackets (see the Supplementary Note 4).
OpenSCAD 1 scripts were used to generate STL meshes for the fabricated structures in this study.

Supplementary Note 2: Simulation
Our framework builds on physical simulation to predict the deformed configuration of the system at a given instant in time. We consider a discrete mechanical model that determines the mechanical behavior of a deformable object based on an elastoplastic model with elastic energy potential W (x,x(t), κ(t)) ∈ R and plastic energy dissipation formulated through the rest configuration update. Here, x(t) ∈ R n is a vector containing n generalized coordinates that spatially discretize the kinematic state of the shell in different configurations,x(t) refers to the undeformed configuration, κ(t) is a vector grouping all material stiffness parameters, and t is the time instant.
Temporal effects are modeled through the explicit dependencies of the undeformed configuration and material stiffness on time. In the next section, we describe how these dependencies are estimated from empirical data. Because our structures morph at low strain rates, we neglect all dynamic effects and take discrete time increments of constant duration, δ = 0.5 s. After each time increment, the undeformed configuration and stiffness parameters are updated quasistatically. We denote the magnitudes corresponding to the k-th time step of the morphing process as x k ,x k , κ k . For notation simplicity, we will drop the superscript corresponding to the morphing time step unless specified. Simulating the static behavior of this mechanical system at the k-th morphing time step implies solving the nonlinear system of differential equations defined by net force equilibrium, i.g., F(x k ,x k , κ k ) = −∇ x W = 0, using standard numerical optimization methods.
Our computational model couples a FEM simulation of the membrane, a rigid body model for bases, and a data-driven spring model for the brackets. In the following sections, we describe the kinematics and mechanics of each of these subsystems separately and then specify how we model the coupling between them and solve the numerical problem.

Discrete kinematics
The bases are modeled using N r prismatic rigid bodies. Each rigid body represents two bases attached to the membrane from the opposite sides. The kinematic state of the i-th rigid body can be determined by the position of the center-of-mass v i ∈ R 3 , together with its rotation r i ∈ R 3 , expressed in angle-axis format. The corresponding rotation matrix R i can be easily computed using the well-known Rodrigues formula where θ = ||r||, u = r/||r|| and [u] × is the cross product matrix of u, i.e., the matrix such that [u] × x = u×x. This allows us to express the position of any point p j in the local coordinates of the i-th rigid body through the non-linear relation All points lying on the surface of a base are coupled to it implicitly. The positions of these points can be expressed in terms of the center-of-mass and rotation of the rigid bodies following the non-linear expression in the Supplementary Equation 1. Therefore, the geometric configuration of the shell can be completely determined by the positions and orientations of the rigid bodies together with the membrane vertices that are not coupled to the rigid bodies. This leads to a total of N t = 6N r + 3N f degrees of freedom which we group in the state vector

Discrete energies
Given this discretization, the mechanical behavior of the shell can be described by a conservative elastic potential W aggregating the contributions of the data-driven springs W c , shear-resisting energy W s , and the membrane W m : here, U ij refers to the unit cell joining the i-th and the j-th bases, T i refers to the i-th element of the membrane discretization, and G i is the deformation gradient of the membrane evaluated at this element. Let us separately explain each of the energy terms: • The data-driven spring energy W c , has the following expression: where L c is the collision distance for the spring determined by the bumpers geometry. Initially, the elastoplastic behavior of the spring follows a data-driven model W d (b, h, t, x) for a given bracket length, b, bracket thickness, h, and time spent under water, t (see next section). It depends on spring deformation, x =L − L(s q ij ), whereL is the rest length of the spring. To account for plastic effects, the rest length of the spring is updated after each time increment following the schemē Once the current length of the spring is smaller than the collision distance, the actuation is stopped and the collision distance is enforced using a soft constraint defined through a high stiffness constant, κ c .
• The shear-resisting energy W s , has the following expression: where L(s ij ) = s i − s j , is the distance between the spring segment end points, R 0 is the difference between distances in the initial morphing time step, and κ s is a constant.
• For the membrane energy W m , we use a classical FEM formulation with an incompressible Neo-Hookean material [4]. Continuum magnitudes are interpolated from nodal values using linear basis functions which allows us to discretely approximate the deformation gradient G = ∇mm. Here, the undeformed configuration can be computed from the membrane state at the initial configurationm = τ −1 m 0 , where τ is the pre-stretch factor.

Coupling and solver
At each time step of the morphing process, we formulate and solve the nonlinear system of differential equations defined by the net force equilibrium −∇ x W = 0, by minimizing the discrete elastic potential in the Supplementary Equation 2. We solve this problem using Newton-Raphson method with Strong Wolfe convergence conditions for step length selection. Solving this problem efficiently requires analytically computing both the first ∇ x W and second ∇ 2 xx W derivatives of the elastic potential. As introduced above, points lying on the surface of the bases, p, are implicitly coupled to the rigid bodies through the Supplementary Equation 1. Hence, these derivatives can be easily computed using the chain-rule: where ∇ x f is a selection matrix of the free membrane vertices and ∇ x p and ∇ 2 xx p can be computed from the Supplementary Equation 1. One common technique in rigid-body simulation to simplify this computation is to keep rotational degrees of freedom r i close to zero. This is done by updating the local coordinates of the coupled points in the base frame after each successful iteration, i.e., for the j-th point attached to the i-th rigid body, p 0 ji ← R(r i )p 0 ji , r i ← 0. We provide the following data as an example of the simulation scales. Our most complex model with self-interweaving shape (Fig. 3d) contains 549 rigid bodies and 7942 membrane elements. The full morphing process simulation (240 time increments) takes 43 minutes in total.

Supplementary Note 3: Material measurement and modeling
As discussed in the Supplementary Note 2, we represent brackets in simulations by data-driven springs and shear-resisting elements. Here we describe our approach to the mechanical modeling of the data-driven components. The brackets in our structure undergo large deformations and are made of a material with nonlinear elastic properties and time-dependent softening. This combination of material and geometric nonlinearities leads us to a data-driven effective spring model. We performed all measurements in settings that resemble conditions brackets are subjected to in an assembled structure. We first formulate the datadriven elastoplastic spring model W d (b, h, t, x) (introduced in the previous section) and then describe our fitting strategy. We discuss several polynomial fittings in this section and display their output units in square brackets.

Physical model of brackets
Our elastoplastic bracket model is motivated by a set of material tests described below. The elastic component is expressed through an effective stiffness, and the plastic behavior is described as a dissipation of internal elastic energy due to deformation. We make the approximation of assuming a constant plasticity fraction, η, which we obtained experimentally. The plastic part of the displacement is then given by x pl = ηx and the elastic part is x el = (1 − η)x.
We aim to obtain elastic energy formulations W d (b, h, t, x el ) that are functions of time t spent in hot water for each valid combination of bracket length b and thickness h. These formulations are modelled as trilinear interpolations between polynomials p b,h,t (x el ) that are defined on a regular grid.
Under the assumption of a constant plasticity fraction and with material properties corresponding to time t, we relate elastic energy W d (b, h, t, x el ) to the total external workW done through monotonic bracket displacements x in the following manner: Here,W are trilinear interpolations of polynomials p work b,h,t (x) [Nm] that are similarly related to p b,h,t (x el ): These polynomials are obtained by integrating force over displacement: where p load b,h,t (x) [N] are polynomials representing the loads exerted on bracket over applied displacements. The polynomials p load b,h,t (x) are fourth-order in x with no free term. The fitting methods for obtaining these polynomials are discussed in the following subsection.
It is challenging to measure time-dependent force-displacement relationships directly in an experimental setup since specimen submersion and loading takes a significant amount of time relative to our actuation time ranges. Specimen submersion and loading takes 8 seconds while the target deformations last approximately 30 to 80 seconds. This restrains us from an assumption that material properties are "fixed" at time t for the material measurements. Additionally, we collected a higher density of data in time rather than in displacement due to the capabilities of our experimental setup. This leads us to dividing the fitting problem into simpler components. We first fit (inverted) displacement-force relationships and then use the obtained model to reconstruct the desired force-displacement model p load b,h,t (x). Our displacement-force model is the following: where p dry [m] is a third-order polynomial (limited to first order in b and h) for which all terms have F as a multiplier and p wet [log mN −1 s −1 ] is a fourth-order polynomial. Polynomial p dry represents displacementforce relationships before putting into water while p wet is the logarithmic evolution of deformation rates, divided by load, in water. The latter formulation restricts the deformation speed to be always positive which is a desirable property for our model and assures zero displacement under zero load. The integral does not generally have a simple explicit representation. We apply the transformation mentioned above to the measured deformation speeds and fit the resulting data The logarithmic function improves the quality of fitting by reducing extreme variation in deformation rates. Note that the formulation of displacement-force relationships is defined as a single continuous function across all parameters in our setup. This allows us to build a consistent model of data collected from all experimental measurements.  Figure 4b) to enforce uniaxial movement. Three types of experiments were conducted: compression tests in both dry and wet states, as well as effective plasticity measurements. For the dry compression tests, displacement-force relations for brackets that hadn't been immersed in water were obtained quasistatically. For the wet compression tests, we immersed specimens into hot water and immediately applied constant loads. These specimens deformed gradually over time since exposure to hot water causes them to soften. For the effective plasticity tests, we compressed brackets by a prescribed displacement, unloaded them, and measured the restoration to derive the plastic component of the deformation. For all tests we used the same sampling of central lengths l (from the

Fitting polynomials with constrained derivatives
Empirical knowledge (such as the fact that thicker brackets deform at slower rates) was used to derive constraints on fitted functions. Similar to standard polynomial regression, we solve a quadratic programming problem y = Xβ + , but instead of β = (X T X) −1 X T y we solve: where A expresses row-wise derivative constraints per point. By providing sufficiently many points in A we significantly improve fitted curve quality. We impose these constraints at all given data points. All the constraints used for each polynomial fitting are listed below.

Compressive loading of dry specimens
This experiment is conducted once for each specimen by increasing the applied load quasi-statically from 0 to 10 N. We fit to this data a polynomial p dry expressing resulting displacements x given the loading F . Note that displacement is a monotonic function of bracket length b, unlike strain. We impose the following derivative constraints: The resulting curves are in the Supplementary Figure 5.

Compressive loading of specimens in water
The next step is to model the evolution of the displacement-force curves over time spent in water. In order to do so, we test the behavior of each specimen under various loads (which are held constant throughout an individual test).
Since the initial force-displacement relationships p dry are known with relatively high precision, we use those curves as our model for t = 0. We fit the displacement rate data to reconstruct the time-evolution of bracket behavior in water. We impose the following constraints on the derivatives:  Figure 6. We filter out displacements larger than 2 /3 of the initial length since brackets may be damaged and display inconsistent behavior at that point.
After obtaining a displacement-force relationship, we invert it with respect to displacement. We densely resample our target parameter space of b, h, t and refit force-displacement curves represented by polynomials p load b,h,t (x) with constrained derivatives: Effective plasticity tests In order to properly simulate actuation of the structure, we examine elastic energy dissipation during bracket compression. We approximate this phenomenon through an elasto-plastic model. We do not study the dependence of plasticity on temperature because all experiments occur at 56 o C. Thus, we analyze the dependence of plasticity on deformation rates. This experiment is done for specimens with fixed chosen parameters l = 8 mm and h = 0.45 mm. We do multiple loading tests with different deformation rates and measure the restored strain after unloading. We observed that there is no significant dependence on deformation rates (see the Supplementary Figure 7). We did not observe consistent strain restoration across different bracket thicknesses and lengths, with mean plastic fractions in a range 15%-25%. We use a constant mean plasticity fraction η = 20% in our simulations.

Supplementary Note 4: Temporal programming
In our design pipeline, the user specifies a time landscape. This is a smooth scalar field over the target surface that represents the desired deformation completion time at each point. In our implementation we define it as a piecewise linear function on top of the flat stencil triangulation. For each unit cell, we compute the average desired actuation time, t * , using values specified at the ends of the associated stencil edge. Then we use t * to configure the thicknesses of the brackets associated with the unit cell. This computation is based on the effective force-displacement curves modelled by trilinear interpolation of polynomials p load b,h,t (x), which we denote as F (b, h, t, x). Each unit cell's deformation is consistent with its neighbors' so long as the time landscape is sufficiently smooth.
Both pairs of brackets on the opposite sides of the unit cell have the same initial length b 0 by construction. We first compute their length in the fully actuated state b 1 (in general, different for different sides). Then, for the target actuation time t * , we can formulate the bracket thickness h configuration problem as follows: where F mem is an approximation of the traction generated by the membrane dependent on the initial unit cell length, b 0 , and on a parameter describing membrane deformation, d 1 with its value at the actuated state. Here we use the force balance between the brackets and the membrane at the target deformation after time t * spent in water, ignoring plasticity (Supplementary Figure 8b).
We build a linearized model of the membrane, representing it as a segment spring connecting the centroids of the end points of all brackets in a unit cell, and setting F mem (b, d) as a first-order polynomial in each of the variables. Here b mimics an initial membrane spring length and d is its displacement. We fit this model by sampling initial bracket lengths and deformations. It captures membrane forces well in our setting (Supplementary Figure 8a).
Our approach to finding the thickness h is a binary search through our interpolated model of F (b 0 , h, t * , b 0 − b 1 ) to match a known value of F mem (b, d).
We need a specific treatment of boundary unit cells since the membrane ends there and has less stretch initially. To account for this, we reduce the force evaluation by a factor of 0.6 which was empirically found by comparing simulations of our examples with the membrane represented by linear springs to FEM membrane simulations.
Once the desired time landscape is specified and bracket thicknesses are computed accordingly, we start the simulation process described in the Supplementary Note 2. Unit cells' actuation completion might deviate from the specified time landscape due to imperfections in the linearized membrane force estimation and plastic effects. Apart from that, the resulting morphing process might not reach the final goal, for example there might be collisions on the way. In these cases time landscape has to be edited and the simulation of morphing has to be recomputed until the goal is achieved.

Supplementary Note 5: Fabrication procedure
We developed a custom fabrication procedure for the shells. It is described in detail in this section.
1. 3D print the outer layers of the shells using a Stratasys J750 printer. We also printed an outer frame to reduce undesired deformation of the structure during fabrication. Printing time for all of our shells is approximately 2 hours (it mostly depends on the number of layers, which is fixed). It is twice as long for the first petalled shape (Fig. 1c) since both sides do not fit on a single printing tray. We remove all overhangs, which does not affect shape replication but makes the cleaning process drastically easier: we only quickly remove almost all support material by several shaving movements with a sharp scraper and airflow the lattice to erase the rest.
Supplementary Figure 10: Membrane stress relaxation over the course of 24 hours. Evolution of the force generated by a dog-bone membrane specimen under a constant stretch factor of 3.
2. Use a star-shaped device to manually stretch a latex sheet uniformly by gradually wrapping it around and fixing at the device's teeth (Supplementary Figure 9a). Friction between the membrane and the device is sufficient to keep the membrane stretched without any additional clamps. A stretch factor of 3 (900% the area) is enforced by matching markers on the sheet to the tips of the device's teeth. Note that large deformations of the membrane lead to slight stress relaxation in the long run. Given the variability of the rest of the system and our simplified membrane model, we neglect this effect. It is however possible to take it into account since after approximately 1.5 hours stress relaxation does not progress much further (Supplementary Figure 10).
3. Clean both surfaces of the sheet with 2-Propanol for better gluing. Apply super glue on a plastic film and distribute it uniformly with a brush. Transfer glue from the film to the 3D printed structures through contact with the plastic film (Supplementary Figure 9b). 4. Glue one lattice to the membrane and pass push pins through the alignment holes on the lattice and through the latex membrane (Supplementary Figure 9c). The membrane does not rip since we pinch it through an isolated area.
5. Glue the opposite side to the latex sheet matching the holes with corresponding alignment pins.
6. Distribute additional super glue on the latex membrane in the region surrounding the lattice perimeter (Supplementary Figure 9d).
7. Wait 5 minutes and then release the latex sheet to a stretch factor near 2. Then using a scalpel cut out the structure, keeping a tiny amount of degenerated latex at the border to prevent membrane ripping (Supplementary Figure 9e).
8. Submerge the structure in 56 • C water and wait for morphing to complete.
9. Take out the deformed shape and let it dry. Under normal conditions the drying process may take roughly 15 minutes.
The whole fabrication process after 3D printing and before drying takes 30-50 minutes.

Supplementary Note 6: Mechanical measurements of shells
We performed a set of mechanical measurements of the shells for estimation of their load-bearing capabilities. We fabricated rectangular and flat regularly tessellated shells composed of identical unit cells of length 7 mm, bracket thickness 0.5 mm, and bumpers 0.3 mm with total dimensions 60 × 55 mm. Since the shell structure is not isotropic, we performed the tests aligning the deformation to one of the three axes parallel to the edges of the hexagonal pattern. The resulting plots are shown in the Supplementary Figure 11. In all tests the deformation speed was approximately 1 mm/min. It is intuitive that the shells have higher resistance to pure in-plane compression due to the bumper contacts in contrast to the other cases when only the membrane and the brackets are loaded.