Abstract

Smart materials provide a means by which we can create engineered mechanisms that artificially mimic the adaptability, flexibility, and responsiveness found in biological systems. Previous studies have developed material-based actuators that could produce targeted shape changes. Here, we extend this capability by introducing a novel computational and experimental method for design and synthesis of a material-based mechanism capable of achieving complex pre-programmed motion. By combining active and passive materials, the algorithm can encode the desired movement into the material distribution of the mechanism. We use multimaterial, multiphysics topology optimization to design a set of kinematic elements that exhibit basic bending and torsional deflection modes. We then use a genetic algorithm to optimally arrange these elements into a sequence that produces the desired motion. We also use experimental measurements to accurately characterize the angular deflection of the 3D-printed kinematic elements in response to thermomechanical loading. We demonstrate this new capability by de novo design of a 3D-printed self-tying knot. This method advances a new paradigm in mechanism design that could enable a new generation of material-driven machines that are lightweight, adaptable, robust to damage, and easily manufacturable by 3D printing.

Introduction

Robotic mechanisms are among the most challenging systems to engineer due to their highly complex, multidisciplinary, and dynamic nature. However, biological organisms have mastered this task using a variety of approaches. This is one of many domains in which nature serves as a source of inspiration and a benchmark against which we measure our progress. Historically, engineered robots have contained rigid links connected via complex actuators such as electric motors [1]. This combination of rigidity and component complexity makes these systems susceptible to failure and also makes them less maneuverable, particularly in tight spaces. By contrast, analogous mechanisms found in nature, such as an elephant trunk or an octopus tentacle, exhibit fluid motion and can be contorted to attain a much wider variety of possible configurations (Fig. 1).

Fig. 1
(a) The Canadarm2 robotic arm positioning an astronaut in space; (b) two elephants contorting their trunks; and (c) an octopus with its tentacles in a complex configuration. (All photos are in the public domain)
Fig. 1
(a) The Canadarm2 robotic arm positioning an astronaut in space; (b) two elephants contorting their trunks; and (c) an octopus with its tentacles in a complex configuration. (All photos are in the public domain)
Close modal

By replacing traditional rigid links and actuators with smart materials, we can ultimately create mechanisms whose robustness and maneuverability more closely resemble that of natural organisms. Smart materials exhibit shape changes in response to changes in their environment, such as temperature changes. By exploiting this capability, our study seeks to advance a new paradigm in mechanism design. The resulting mechanisms will have no electromechanical components and no central processing station. Therefore, they can be rapidly and cheaply manufactured via 3D printing. This form of 3D printing in which the fabricated object is designed to change shape over time is sometimes referred to as 4D printing [2]. Additionally, because the motion is no longer restricted to a small number of discrete actuators, 4D-printed mechanisms can have a high concentration of degrees-of-freedom, similar to elephant trunks, which can have 40,000 muscles with no bones or joints, thus enabling highly complex motion. This allows for increased maneuverability and much more complex motion.

Recently the concept of physical artificial intelligence has been introduced as a counterpart to digital artificial intelligence [3]. Physical AI systems rely heavily on the use of smart composite materials to create next-generation robots that are akin to biological organisms. While researchers have made great strides in digital AI for robotics over the past few decades, advances in the development of robots’ bodies, materials, and morphology have not kept pace [3]. However, recent advances in 4D printing and computational design suggest that this emerging research area will have a significant role to play in filling this void.

Several early studies on 4D printing combined multiple smart materials, typically shape-memory polymers (SMPs), to create hinge mechanisms that provided in-plane bending actuation. Ge et al. employed this strategy to create 4D-printed active origami composites [4]. They were able to combine multiple hinge components to create mechanisms such as a self-assembling box and a mock paper airplane, both of which emerged from the 3D printer as a flat sheet before actively morphing into their respective 3D shapes. Raviv et al. used a similar approach in which they systematically assembled 4D-printed actuators, which they referred to as primitives, to create long chains that could change shape and spell out pre-programmed acronyms [5]. Here, the researchers used an algorithm to determine the sequence of primitives necessary to obtain the desired end shape. More recently, Gu et al. used design optimization to create a soft prosthetic hand containing fiber-reinforced elastomeric composites [6]. Unlike the 4D-printed examples cited above, motion at the hand’s joints was achieved via pneumatic actuation. By replacing the electric motors found in conventional prosthetics, the prosthetic hand achieved significant cost and weight savings. Each of these examples demonstrates the power of material-driven actuation, but they maintain key features of traditional robotic mechanisms, namely the use of rigid struts separated by a small number of deformable joints. Consequently, the number of achievable configurations for each designed mechanism remains limited. A more recent study by Gladman et al. took a biomimetic approach to create 4D-printed plant-inspired designs that morph into target shapes encoded into their material architectures using composite hydrogels with controlled anisotropic swelling [7].

Topology optimization provides a computational framework for the design of smart mechanisms with complex programmable motion. This technique was originally developed in the late 80s as a means of generating maximally stiff load-bearing structures [8]. Since its introduction, the method has been adapted to create a wide range of mechanical systems including compliant mechanisms [9], aeroelastic structures [10], cellular materials [11], and active composites [12]. Topology-optimized smart composites contain a combination of active and passive materials distributed freely throughout the volume of the mechanism [13], such that there is no distinction between joints and structural members. Hence, the entire body of the mechanism acts as an actuator, thereby allowing for complex free-form motion, similar to that which we see in nature. Several studies have proposed the use of topology optimization to automate the process of distributing active and passive material phases at the scale of the voxel, with 3D printing being used to fabricate the highly intricate designs that result [12,14,17]. In this way, the movement of the mechanism is programmed into the material distribution. However, the computationally intense nature of the topology optimization method has limited the complexity and programmability of the active mechanisms designed in these earlier studies. In each of the previously published examples, the motion of the mechanism was characterized by small shape changes [14] or rotation about a single axis (i.e., pure bending [13,14] or pure torsion [12,13]).

Here, we introduce a novel method that incorporates topology optimization into a broader hierarchical framework in which the algorithm selects the material layout at two scales: at the smaller scale the material is selected voxel by voxel, while at the larger scale, the algorithm selects the arrangement of a series of optimized kinematic units. This approach allows us to generate large, complex motions while keeping the computational cost tractable. We extend the capability achieved in earlier studies by enabling complex motion characterized by large nonlinear deflections that include rotations about multiple reference axes. Thus, our designs combine the robustness and manufacturability of earlier 3D-printed material-driven mechanisms with the programmability of conventional robotic mechanisms. Additionally, with the entire body of the mechanism providing sensing, actuation, and structural rigidity, this design method represents an important step toward creating engineered mechanisms that can compete with natural organisms in terms of adaptability and free-form morphability.

Methods

The proposed design method proceeds in three stages: (1) topology optimization of basic kinematic elements, (2) experimental characterization, and (3) computational assembly of the mechanism (see Fig. 2). In the first stage, we use multimaterial topology optimization [16] to design a series of kinematic elements that produce fundamental displacement outputs such as bending and twisting. The overall mechanism is an assembly of unitary elements that comprise three distinct classes: bending elements, torsional (i.e., twisting) elements, and neutral elements, whose shape remains fixed during the activation of the mechanism. Together, these elements enable rotational displacement about all three coordinate axes in both the positive and negative directions. Prior to actuation, each element is in the shape of a rectangular prism, and the optimization algorithm is used to distribute two smart materials throughout this prescribed volume in order to maximize the element’s angular displacement. By dividing the design task into two stages (i.e., topology optimization and mechanism assembly), we significantly reduce the computational cost of the problem. Performing time-dependent three-dimensional finite element analysis (FEA) and optimization of the full mechanism in a single step would be highly computationally expensive. Our approach reduces the number of degrees-of-freedom in the FEA model by more than an order of magnitude.

Fig. 2
The hierarchical design framework: (a) structural topology optimization of kinematic elements; (b) the thermomechanical programming cycle used to trigger the shape-memory effect in the polymeric design materials. The specimen is initially heated to a hot temperature TH, which is above the glass-transition temperatures of both design materials (Tg1 and Tg2). It is then stretched and cooled to a temperature Tc, which is below both glass-transition temperatures. It then undergoes a relaxation period during which the stretching force is removed. Lastly, it is reheated to a temperature between Tg1 and Tg2 to trigger the desired motion. Note that the blue and yellow images represent the topology-optimized bending beam, with the dashed blue outline indicating the original shape and size prior to thermomechanical loading. The displacement is then measured via digital image correlation. (c) In the final stage of the process, we perform a second optimization procedure to obtain the required sequence of elements before fabricating the full mechanism.
Fig. 2
The hierarchical design framework: (a) structural topology optimization of kinematic elements; (b) the thermomechanical programming cycle used to trigger the shape-memory effect in the polymeric design materials. The specimen is initially heated to a hot temperature TH, which is above the glass-transition temperatures of both design materials (Tg1 and Tg2). It is then stretched and cooled to a temperature Tc, which is below both glass-transition temperatures. It then undergoes a relaxation period during which the stretching force is removed. Lastly, it is reheated to a temperature between Tg1 and Tg2 to trigger the desired motion. Note that the blue and yellow images represent the topology-optimized bending beam, with the dashed blue outline indicating the original shape and size prior to thermomechanical loading. The displacement is then measured via digital image correlation. (c) In the final stage of the process, we perform a second optimization procedure to obtain the required sequence of elements before fabricating the full mechanism.
Close modal

The mechanism is composed of two shape-memory polymer materials. Although both materials contain shape-memory properties, one of the design materials is chosen to have a higher transition temperature so that when the mechanism is heated to activate the shape-memory effect, this higher-temperature material remains unactivated. The resulting difference in the displacement response within the activated and unactivated material regions is harnessed to produce intricate three-dimensional motion. Therefore, in this context, we refer to the material with the lower transition temperature as the active material, and we refer to the material with the higher transition temperature as the passive material. Within the optimization algorithm, we use finite element analysis to simulate and predict the transient displacement response of each kinematic element. This information is then passed to an optimization algorithm to update the material distribution within the mechanism. This process is carried out iteratively until it converges to an optimal design. In this way, the algorithm systematically determines which material (active or passive) should be used to populate each voxel within the overall volume. The output of the algorithm is a voxel-by-voxel description of the desired material distribution. A detailed description of the topology optimization procedure is provided later in this section.

In the second stage of the process, we experimentally measure the displacement response of each kinematic element. For this task, we 3D print a representative prototype of each element class and measure the angular displacement of the element in response to the designed thermomechanical cycle. Although the optimization algorithm uses numerical simulation, which offers an approximate prediction of the displacement response, this simulation can be subject to modeling errors. By performing this measurement step, we can accurately characterize the anticipated displacement response of each kinematic element. In this way, we effectively negate any errors associated with the computational model, since only the experimental measurements are passed to the final stage of the process, where a second algorithm assembles the elements into a specific sequence.

To trigger the displacement response, each element must undergo a three-step thermomechanical programming cycle. In the final step of the cycle, we heat the element to an elevated temperature to activate its shape-memory effect and produce a self-morphing motion. The resulting displacement is then measured using a digital image correlation procedure. For this measurement, we imported the image file into the matlab software environment. From here, we manually selected key locations on the image, and we used a matlab script to extract the xy-coordinates of these points and compute the relevant displacement angles. Note that although we use a liquid bath for heating and activating the mechanism, this approach may not be practical for certain applications. The general methodology presented here is also compatible with other forms of heating, such as joule heating, which has been successfully applied to 3D-printed polymers that have undergone a carbonization process to increase their electrical conductivity [17].

In the final stage of the design process, we use a genetic algorithm to determine the sequence of kinematic elements required to achieve the target displacement. In this study, we have chosen to design a self-tying knot as a demonstration example. This problem is particularly challenging because the knotted path combines bending and torsion about multiple reference axes. The knot also exhibits large geometrically nonlinear displacements, which require a fully three-dimensional kinematic model, similar to that which is used in robotics problems [18]. We must also implement design constraints to prevent the various sections of the knot from colliding with one another during the motion. The genetic algorithm starts out by generating random sequences of kinematic elements and calculating the displacement trajectory of each sequence using a forward kinematics model. Each sequence serves as a candidate design for the self-tying knot mechanism. The calculated final shape of the sequence is then compared with the shape of the ideal knot, to quantify the fitness of each design (see Fig. 3). After successive iterations, the algorithm converges to produce a sequence whose end shape is as close as possible to that of the ideal knot. Below we provide a detailed description of the computational and experimental procedures implemented during this study.

Fig. 3
The designed self-tying knot (red) superposed onto the ideal knot (green). To find the optimal sequence of kinematic elements, we select two anchor points distributed along the length of the knot. Each point has an analogous point on the ideal knot located at a prescribed arc length from the root of the knot. The algorithm searches for the element sequence that causes the position and orientation of the designed knot to be as close as possible to the position and orientation of the ideal knot at the designated anchor points.
Fig. 3
The designed self-tying knot (red) superposed onto the ideal knot (green). To find the optimal sequence of kinematic elements, we select two anchor points distributed along the length of the knot. Each point has an analogous point on the ideal knot located at a prescribed arc length from the root of the knot. The algorithm searches for the element sequence that causes the position and orientation of the designed knot to be as close as possible to the position and orientation of the ideal knot at the designated anchor points.
Close modal

Topology Optimization.

To obtain the optimal material distribution within the kinematic elements, we perform a multimaterial topology optimization. Within the optimization algorithm, the thermomechanical response of the material is simulated using the FEA. The FEA model for the 3D torsional unit contains a mesh of 625 eight-node hexahedral elements. The FEA model for the bending element used a 2D mesh containing 2880 4-node quadrilateral elements. The optimal design was then extruded to achieve a cross section with an aspect ratio of 1. Both models assume geometrically linear, small deformations. The material model is based on the small-strain shape-memory polymer constitutive model proposed by Baghani et al. [19]. To solve the numerical optimization problem, we use the gradient-based method of moving asymptotes (MMA) [20]. The design sensitivities required by the MMA algorithm are computed analytically using a transient adjoint formulation. Further details of the FEA model and the sensitivity analysis procedure can be found in Ref. [13].

There are two classes of kinematic elements that must be designed via topology optimization: the bending element and the twisting element. Once we have the material distribution for each class of element, we can simply rotate or flip the design to obtain angular actuations for all three rotational degrees-of-freedom, as shown in Fig. 4. Note that during the final synthesis phase of the design framework, the algorithm will select from among six actuators at each point along the chain. However, all four bending actuators will have the same material distribution, and both twisting actuators will have the same material distribution. In addition to these six options, the optimizer may also insert a neutral element, whose shape remains unchanged during the activation process. This element can be used as an extender to shift the mechanism's kinematics along a given axis, and it is included in Fig. 4(a) to illustrate the initial shape of all actuators prior to actuation. The topology optimization procedure is used to obtain a precise distribution of the active and passive materials that will produce the desired rotational motion at the end of the thermomechanical programming cycle. For this task, the authors have developed a novel material representation scheme, which is an extension of the Solid Isotrpic Material with Penalization (SIMP) method, which is commonly used for elasticity problems [21]. In our case, we must continuously interpolate the elastic modulus, E, but also the viscosity coefficients, ηr and ηg (one for the glassy and rubbery phases of each design material), and the thermal expansion coefficients, α. Hence, within each element, the effective value of a given material property, Ψ, is computed as an interpolation between the actual properties of the two active and passive design materials, Ψ1 and Ψ2, respectively. The interpolation formula is given in Eq. (1)
(1)
where ρ is the design variable assigned to the element in the FEA model. When ρ = 0, the element will have the properties of material 1 (i.e., Ψeff = Ψ1), and when ρ = 1, the element will have the properties of material 2. The parameter p is a constant that is used to penalize hybrid material states. When p = 3, this pushes the optimization search toward binary designs. Under our interpolated scheme, penalization is applied only when interpolating the elastic moduli, E. When interpolating all other material properties, the penalization is turned off by setting p = 1. Note that the use of Eq. (1) for material interpolation is an extension of the commonly used two-phase SIMP interpolation scheme [21].
Fig. 4
Actuator options for the kinematic elements, which form the building blocks of the self-tying knot mechanism. (a) Neutral element, (b) rotation about the z-axis in the positive direction (bending element), (c) rotation about the y-axis in the positive direction (bending element), (d) Rotation about the x-axis in the positive direction (twisting element), (e) rotation about the Z-axis in the negative direction (bending element), (f) rotation about the Y-axis in the negative direction (bending element), and (g) rotation about the x-axis in the negative direction (twisting element).
Fig. 4
Actuator options for the kinematic elements, which form the building blocks of the self-tying knot mechanism. (a) Neutral element, (b) rotation about the z-axis in the positive direction (bending element), (c) rotation about the y-axis in the positive direction (bending element), (d) Rotation about the x-axis in the positive direction (twisting element), (e) rotation about the Z-axis in the negative direction (bending element), (f) rotation about the Y-axis in the negative direction (bending element), and (g) rotation about the x-axis in the negative direction (twisting element).
Close modal

For the bending element, all motion occurs within xy-plane. Therefore, we can treat the material distribution problem as a two-dimensional problem and extrude the solution along the z-axis to obtain the full three-dimensional material distribution. Figure 5 shows the geometry and boundary conditions for the two-dimensional design domain used for the bending elements. This figure also shows the interpolated and extruded material distribution, which forms the basis of the 3D computer-aided design (CAD) model that was ultimately fabricated using 3D printing. During the interpolation process, we begin with the material distribution field, represented by the variable ρ, which is a piecewise constant within each finite element. These data are then interpolated to extract the precise location of the material interface to obtain a high-resolution representation of the internal material distribution. In this context, the material interface is selected as the contour or level set corresponding to ρ = 0.5. Note that due to the nature of the polymer cross-linking across the material interface, we assume ideal bonding at the material interface, and therefore, there is no need to include a special model to account for interfacial effects [4].

Fig. 5
Topology optimization of the bending element: (a) Geometry and boundary conditions of the design domain. The design materials must be distributed within this rectangular domain. Note also that the axial force, F, is applied only during the initial stage of the thermomechanical programming cycle. (b) Optimized material distribution in pixel form. The blue region represents the active material and the yellow region represents the passive material. The color bar indicates the value of the design variable ρ. Each pixel in the image represents a single element in the finite element mesh. (c) The interpolated 3D material distribution after extrusion in the z-axis. Note that the interpolation process removes all intermediate-density material appearing along the interface.
Fig. 5
Topology optimization of the bending element: (a) Geometry and boundary conditions of the design domain. The design materials must be distributed within this rectangular domain. Note also that the axial force, F, is applied only during the initial stage of the thermomechanical programming cycle. (b) Optimized material distribution in pixel form. The blue region represents the active material and the yellow region represents the passive material. The color bar indicates the value of the design variable ρ. Each pixel in the image represents a single element in the finite element mesh. (c) The interpolated 3D material distribution after extrusion in the z-axis. Note that the interpolation process removes all intermediate-density material appearing along the interface.
Close modal

The optimization problem statement for the bending element is given in Eq. (2). Note that we seek a material distribution, represented by the design variable, ρ, that will maximize the vertical displacement UyN at the designated node, N (marked with a red dot in Fig. 5, at the end of the thermomechanical cycle ((t=t*)), while constraining the optimizer to ensure that the total volume of active material VSMP1 does not exceed 70% of the total volume of the element. We impose a limit on the volume of the active material because this material is less stiff than the passive material. Therefore, it is necessary to limit the amount of active material available to the optimizer in order to ensure that all kinematic units have the necessary axial stiffness and that all portions of the mechanism are stretched by a similar amount during the thermomechanical programming cycle. The specific choice of a 70% volume fraction was chosen as a tradeoff between achieving a stiff design and a design that could produce large displacements, since the active material is what drives the morphing behavior

(2)

Figure 6 shows the geometry and boundary conditions for the design of the twisting element. This figure also contains the optimized material distribution generated by the algorithm, along with the interpolated material distribution within the 3D CAD model. Note that due to the high computational cost of the three-dimensional transient thermoelastic finite element analysis and the accompanying sensitivity analysis, we use a relatively coarse mesh for the initial topology optimization; however, after the interpolation process, we obtain a sufficiently smooth material interface.

Fig. 6
Topology optimization of the twisting element: (a) Geometry and boundary conditions of the design domain. (b) Optimized material distribution in pixel form. The black region represents the active material and the white region represents the passive material. (c) The interpolated 3D material distribution.
Fig. 6
Topology optimization of the twisting element: (a) Geometry and boundary conditions of the design domain. (b) Optimized material distribution in pixel form. The black region represents the active material and the white region represents the passive material. (c) The interpolated 3D material distribution.
Close modal
Equation (3) shows the optimization problem statement for the design of the twisting element. Note that we seek to minimize the deflection of node 2 (the lower red dot in Fig. 6(a)) in the z direction, thereby pushing this node inward. At the same time, we constrain the z deflection of node 1 to be a positive number. These two displacements combine to create a twisting motion at the free face of the element.
(3)

The two-dimensional topology optimization of the bending element as well as the 3D topology optimization of the twisting element is implemented in c++, with the PETSc library used for parallelization. Both sets of code are published in a public repository.

Experimental Measurement and Characterization.

The bending and twisting structures are fabricated using the Stratasys Objet260 Connex 3D printer. The printer produces a wide range of digital materials by mixing two base materials: a stiff polymer VeroWhite Plus and a soft rubber-like elastomer TangoBlack Plus [7]. All the digital materials exhibit shape-memory behavior, and each of the digital materials has a unique glass-transition temperature (Tg). We printed a subset of the available digital materials and experimentally determined their glass-transition temperatures. Finally, two digital materials, RGD8530 and FLX9895, were selected as the design materials for the bending and twisting structures. These materials were chosen because they have clearly distinct glass-transition temperatures, allowing us to activate targeted sections of the mechanism using a uniform temperature field.

The bending and twisting structures were built with dimensions of 1 × 0.2 × 0.2 cm so that their aspect ratios of 5:1:1 matched the designs obtained computationally. The scale of these unit structures was decided based on the desired size of the knot mechanism, the size of the 3D printer’s build envelope, and the estimated number of elements within the final sequence. Representative structures with these particular dimensions were fabricated and subjected to the thermomechanical programming cycle, as shown in Fig. 7, to determine their displacement response. The equipment used to characterize the displacement kinematics is described below.

  1. Temperature-controlled water bath: This is used to apply a uniform temperature field to the shape-memory polymer samples over an extended period of time.

  2. Tensile device: A global uniaxial strain is applied to the SMP samples during the deformation step of the thermomechanical programming cycle using an adjustable tensile device. The tensile device connects to two loops attached to either end of the morphing mechanism. The loops are made of VeroWhite, a stiff polymer with a high transition temperature. This material choice ensures that the connecting loop remains relatively rigid during the stretching of the mechanism.

  3. Water tank and water circulator with temperature control: The stretched SMP structures are placed inside a temperature-regulated water tank and maintained at a fixed temperature using a water circulator during the reheating step for their motion activation.

Fig. 7
The thermomechanical programming cycle used for the experimental case studies. TH and TL represent the highest and lowest temperatures of the cycle. The values Tgp and Tga represent the glass-transition temperatures of the active and passive SMP materials, respectively.
Fig. 7
The thermomechanical programming cycle used for the experimental case studies. TH and TL represent the highest and lowest temperatures of the cycle. The values Tgp and Tga represent the glass-transition temperatures of the active and passive SMP materials, respectively.
Close modal

The thermomechanical programming cycle consists of three steps: cooling with deformation (C + D) of the specimen, relaxation (R), and reheating (H) of the specimen to obtain the desired shape change. The procedure is illustrated graphically in Fig. 7. Prior to starting the cycle, we heat the mechanism up to a temperature TH, that is above the glass-transition temperatures of both the active and passive SMP materials so that the entire specimen is in the rubbery phase. We then begin the programming cycle by applying a tensile load to stretch the specimen while cooling it to a temperature TL that is below the glass-transition temperatures of both the active and passive materials. This cooling and stretching step lasts 20 min, during which both materials transition to their glassy phase. Next, during the relaxation step, we release the applied forces while maintaining the specimen at its low temperate TL for an additional 20 min. During this step, the internal stresses in the material reduce the zero; however, the material retains its stretched shape due to the strain fixity property of the shape-memory polymers. In the third and final stage of the cycle, we raise the temperature to 34 °C, which is the glass-transition temperature of the active material. In this way, only the active material transitions back to the rubbery phase and therefore seeks to return to its original undeformed shape. The resulting disparity in strain between the two material regions causes the mechanism to exhibit the targeted shape change.

To measure the deflection angle of the bending element, a chain of 10 bending elements was 3D-printed, and the combined structure was subjected to the full thermomechanical programming cycle. The bending angle for a single element was determined by evaluating the bending angle of the combined structure and dividing it by the number of elements. The bending angle was obtained by analysis of the digital image of the activated structure shown in Fig. 8(b). By taking this averaged measurement, we sought to reduce the impact of outliers and manufacturing defects. Similarly, the rotation angles for individual torsional specimens were calculated by subjecting a sequence of three torsional elements to the full thermomechanical programming cycle and calculating the average rotation angle.

Fig. 8
Experimental determination of the twist and bending angles of the kinematic elements: (a) twist angle illustration, (b) bending angle illustration, (c) activated series torsional elements, and (d) activated series bending elements
Fig. 8
Experimental determination of the twist and bending angles of the kinematic elements: (a) twist angle illustration, (b) bending angle illustration, (c) activated series torsional elements, and (d) activated series bending elements
Close modal

For all kinematic elements, the rotation angle of the activated element is proportional to the strain applied during the thermomechanical programming cycle. For the twisting elements, the test samples were stretched to a strain of 12%, which produced a twist angle of 4deg per element as shown in Fig. 8(a). For the bending elements, eight different samples were tested with strains ranging between 9% and 16%. These results were then interpolated to obtain the expected bending angle for a strain of 13%, which was the applied strain of the self-tying knot during thermomechanical programming. This interpolation produced a bending angle of 25deg for each individual bending element. These results (4deg and 25deg) were then used as the respective activation angles for the twisting and bending elements in the forward kinematic simulation described in the following section. Note that the finite element model predicted angular displacements of 3.5deg and 33deg for the torsional and bending elements respectively. This corresponds to an error of 12.5% for the torsional element and 32.0% for the bending element. The large discrepancy observed in the displacement for the bending beam is primarily due to the assumption of small strains in the finite element model. However, this error is mitigated by the fact that only the experimentally obtained displacements are used in the final kinematics model.

It should be noted that the experimental measurements also contain some inevitable errors and variability. In order to maintain the repeatability of the measurements, we must maintain a consistent protocol in the temperatures and durations of the various stages of the thermomechanical programming cycle. We must also control the process parameters of the 3D printing process. This includes things like ensuring that the intensity of the ultraviolet light used to cure the polymers is kept constant across print jobs in order to maintain a consistent cross-linking density.

Mechanism Synthesis.

The displacement and orientation of the tip of the knot are evaluated computationally using a forward kinematics model. The rotation angle of the bending and twisting elements used in the model is based on the experimentally measured rotation angles discussed earlier. An ideal knot shape is used as a target design to guide the algorithm. The path of the ideal knot is given by the following parametric equation
(4)

Note that the parameter t ranges from −3π/4 to 3π/4. The values of the coefficients in Eq. 4 were selected so that the target knot size would remain within the volume of the 3D printer’s maximum build envelope. The forward kinematics model is combined with a genetic algorithm to arrange the twisting and bending elements in an optimal sequence so that the end shape of the trial knot resembles the shape of the ideal knot as closely as possible. When joining the kinematic elements end to end, we assume that the end faces remain planar and orthogonal to the axis of the element. This assumption is justified by the relatively small angular deflections observed in the experiments and is supported by beam theory. Consequently, we assume the elements can be joined seamlessly with no residual stresses caused by incompatible strains at the interface.

Note that because the knot is symmetric, we design only the half-knot and reflect this solution about the root (base) of the knot to obtain the full knot path. The objective function for the optimizer consists of two terms as follows:
(5)
where Perror represents the error between the Cartesian coordinates {x, y, z} of the trial knot and the ideal knot at designated anchor points. In this study, the anchor points are chosen as the midpoint of the half-knot and the tip of the half-knot. The second term in Eq. 5, Qerror, represents the error in the orientation of the two knots at the tip.
(6)
The coefficients were set to C0 = 1.0 and C1 = 5.0 in order to appropriately scale the two error terms. The subscripts m and t refer to the midpoint and the tip of the knot, respectively. The subscript 0 refers to the coordinates of the ideal knot and the constant coefficient was set to wm = 5.0 to assign a greater weight to the midpoint error, since errors at the midpoint will magnify errors at the tip of the knot. The angles ϕ and θ capture the orientation of the knot segment expressed in spherical coordinates. Note that for the 13-segment half-knot (12 elements plus one neutral element), the midpoint of the ideal knot was approximated using {x, y, z} coordinates evaluated at t = (3π/4)/(12/2 + 1). The forward kinematics algorithm loops over all the elements of the trial knot, evaluating the rotation matrices along each axis (Rx, Ry, and Rz) within a local reference frame as follows:
(7)
(8)
(9)
where α, β, and γ are rotational angles undergone by a 3D body about the x-axis, y-axis, and z-axis, respectively. The total rotation (Re) of the 3D element is computed from the rotation along each of the axes. This is then used to evaluate the global rotation of the full sequence (R)
(10)

For each element in the sequence, offset distances along each axis (δx, δy, and δz) are calculated for the displaced kinematic element.

For a bending element, shown in Fig. 9, the offset distances are given as follows:
(11)
where L is the length of the element in the x-dimension prior to thermoelastic deformation. The offset distances for the twisting element are given below.
(12)
Fig. 9
Offset distances shown on a single activated kinematic bending element
Fig. 9
Offset distances shown on a single activated kinematic bending element
Close modal

The offset distances are used to evaluate the global change in position vector (ΔX) which is then used to calculate the position of the midpoint and tip of the full trial knot. The steps followed in the forward kinematics framework are shown in Algorithm 1.

Pseudocode for the forward kinematics analysis

Algorithm 1:

R = eye(3, 3) /* initialize rotation matrix as the identity (i.e., no initial rotation)*/

X = zeros(3, 1) /* initialize the position vector*/

fori ← 1, 2,…., Ndo

/* loop over all the elements of the sequence */

/* Evaluate Rx, Ry and Rz given by Eqs. (7)(9) */

/* Form a vector with offset distances δx, δy and δz using Eqs. (11) and (12) */

ΔX=R(δxδyδz) /* Evaluate change in position vector ΔX */

Re = RzRyRx /* Evaluate the total rotation based on the current element using Eq. (10) */

RRRe /* Update the total rotation matrix R */

XX+ΔX /* Update the position vector X */

The forward kinematics algorithm is combined with matlab’s genetic algorithm function “ga.” The forward kinematics module evaluates the objective function for each trial knot design, which is then fed into the optimizer. The optimization algorithm was initially allotted 10 elements with which to construct the self-tying knot. If the optimization failed to generate a complete knot (where the tips penetrate the opposing loops on either side), then the number of elements was increased by one. In addition to requiring a complete knot, the algorithm also rejected any designs that caused collisions in which two segments of the knot would come into contact with one another at any point during the morphing process. To detect potential collisions the morphing motion was decomposed into 25 increments, and the Euclidean distances between all node pairings were evaluated at each increment.

Results and Discussion

The minimum number of elements with which the algorithm successfully produced a complete self-tying knot design with no collisions was 13. The optimized sequence for the half-knot is given in Fig. 10. In the figure, the identifiers ag correspond to the letters of the seven kinematic element options shown in Fig. 4.

Fig. 10
Offset distances shown on a single activated kinematic bending element
Fig. 10
Offset distances shown on a single activated kinematic bending element
Close modal

As indicated in Fig. 10, the two halves of the designed knot are 3D-printed together, and they share a single connecting loop on both ends, for simultaneous thermomechanical programming. After stretching and relaxation, the connecting loops are removed with scissors, and the two halves of the knot are mounted on a specialized stand prior to being placed in the heated water tank for final activation. The stand contains a two-sided flexible elastomer socket, and the two half-knots are slid into either side of the socket at End 1 as indicated in Fig. 11.

Fig. 11
CAD model of the mounting apparatus used to hold the mechanism in place during activation (note that this image shows only the back half of the stand so that its cross section is visible)
Fig. 11
CAD model of the mounting apparatus used to hold the mechanism in place during activation (note that this image shows only the back half of the stand so that its cross section is visible)
Close modal

The activated self-tying knot shape is compared with the ideal knot shown in Fig. 12. The blue design represents the ideal knot generated mathematically, superimposed on the computationally generated self-tying knot shape, shown in green.

Fig. 12
Comparison of the converged self-tying knot shape (bright) with the mathematically generated ideal knot shape (dark)
Fig. 12
Comparison of the converged self-tying knot shape (bright) with the mathematically generated ideal knot shape (dark)
Close modal

Note that the geometry shown in Fig. 12 does not account for the effects of gravity. While the mechanism remains submerged in the water bath, the gravitational force is counteracted by the buoyant forces applied by the water and therefore we observe minimal sagging of the knot mechanism. Consequently, the photos of the submerged 3D-printed knot (Fig. 14(b)) show good agreement with the forward kinematic simulation represented by Fig. 12. However, once the knot is removed from the water, gravitational effects become more apparent and the knot sags, as shown in Fig. 14(a). Figure 13 shows a computational simulation of the deformation of the activated knot mechanism when subject to gravity with no counteracting buoyant force. The simulation was performed using the commercial software package ANSYS with 13 space frame elements (one for each element in the forward kinematics model). Each element was assumed to be linearly elastic with 12 degrees-of-freedom per element. The shape of the deflected knot mechanism as determined by the simulation shows good agreement with the experimental results shown in Fig. 14(a).

Fig. 13
The activated self-tying knot with elastic deflection due to gravity (red) shown with the original undeflected knot (green)
Fig. 13
The activated self-tying knot with elastic deflection due to gravity (red) shown with the original undeflected knot (green)
Close modal
Fig. 14
The final 3D printed self-tying knot mechanism: (a) the end shape of the self-tying knot shown from different viewing angles; and (b) time lapsed photos of the self-tying knot during the activation phase when the knot is submerged in a heated water bath. The images are separated by three-second intervals, with the full motion occurring over 24 seconds. Note that in (a) the mechanism has been removed from the water and is no longer subject to buoyancy forces, therefore its shape sags due to gravity.
Fig. 14
The final 3D printed self-tying knot mechanism: (a) the end shape of the self-tying knot shown from different viewing angles; and (b) time lapsed photos of the self-tying knot during the activation phase when the knot is submerged in a heated water bath. The images are separated by three-second intervals, with the full motion occurring over 24 seconds. Note that in (a) the mechanism has been removed from the water and is no longer subject to buoyancy forces, therefore its shape sags due to gravity.
Close modal

The disparity between the target shape and the observed shape of the 3D-printed knot has several contributing factors. The image-based measurement of the angular displacements of the 3D-printed kinematic elements has some inherent error, which accumulates as the number of elements in the sequence increases. Additionally, the genetic algorithm used to determine the optimal sequence of elements has a limited number of degrees-of-freedom with which to approximate the target shape. Therefore, it is necessary to select the number of kinematic elements that will offer an acceptable tradeoff between the accuracy of the knot shape predicted by the kinematic model and the precision of the genetic algorithm.

Figure 14 shows the final shape of the actual 3D-printed self-tying knot mechanism, along with a series of time-lapsed photos of the mechanism to illustrate the full range of motion.

Conclusions

The proposed design framework combines several novel features that enable unique functionality. Whereas a single-step approach could potentially lead to a computationally intractable optimization problem due to the high dimensionality of the design space, our hierarchical approach allows for efficient exploration of the vast design space, which provides for enhanced programmability. The resulting mechanisms can therefore generate complex motion, all of which is systematically encoded into the material distribution. This capability is also enabled by our hybrid approach, in which experimental measurement is embedded into the computational design framework. Second, by replacing electromechanical components with material-based actuation, we make the mechanisms more lightweight, more robust to failure, and readily manufacturable via 3D printing. Furthermore, these material-based mechanisms are highly miniaturizable. This could facilitate the creation of micro- and nanorobots, which are particularly useful for in vivo biomedical applications such as targeted drug delivery [22]. Lastly, the lack of a central processor means that both sensing and actuation occur locally, allowing for rapid response to environmental changes. Future work will expand the proposed method beyond chain-like mechanisms to include 2D and 3D arrays of unit structures.

In their 2012 book Fabricated, Kurman and Lipson [23] speculate that 3D printing will revolutionize the way we procure everyday household items such as toothbrushes, which can now be printed at home on demand. Advances in 4D printing may take this a step further by allowing consumers to 3D print previously unseen items such as self-tying shoelaces and combs with adaptive bristle density. Advances like the ones presented in this study could make possible an entirely new class of technologies that rely on programmable material-based robotic systems. More importantly, this capability could provide a key ingredient necessary for realizing aspirational life-saving technologies such as artery-clearing microrobots, self-tying sutures, and many other disruptive technologies that today seem unimaginable.

Footnotes

3

See Note 2.

4

See Note 2.

Acknowledgment

This work was funded by the National Science Foundation (Grant No. 1663566).

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The data and information that support the findings of this article are freely available.2

Materials and Correspondence

Correspondence and inquiries should be addressed to Kai A. James (kai.james@gatech.edu). All computer code used to generate the results presented in this article is publicly available for download at the following repository.3

Statement on Competing Interests

The authors declare no competing interests.

Data Availability

All computer code used to generate the results of this study is publicly available via the following online repository.4

References

1.
Denket
,
M. A.
,
2006
,
Frontiers in Robotics Research
,
Nova Publishers
,
Hauppauge, NY
.
2.
Tibbits
,
S.
,
2014
, “
4D Printing: Multi-Material Shape Change
,”
Archit. Des.
,
84
(
1
), pp.
116
121
.
3.
Miriyev
,
A.
, and
Kovač
,
M.
,
2020
, “
Skills for Physical Artificial Intelligence
,”
Nat. Mach. Intell.
,
2
(
11
), pp.
658
660
.
4.
Ge
,
Q.
,
Dunn
,
C. K.
,
Qi
,
H. J.
, and
Dunn
,
M. L.
,
2014
, “
Active Origami by 4D Printing
,”
Smart Mater. Struct.
,
23
(
9
), p.
094007
.
5.
Raviv
,
D.
,
Zhao
,
W.
,
McKnelly
,
C.
,
Papadopoulou
,
A.
,
Kadambi
,
A.
,
Shi
,
B.
,
Hirsch
,
S.
, et al
,
2014
, “
Active Printed Materials for Complex Self-Evolving Deformations
,”
Sci. Rep.
,
4
(
1
), p.
7
.
6.
Gu
,
G.
,
Zhang
,
N.
,
Xu
,
H.
,
Lin
,
S.
,
Yu
,
Y.
,
Chai
,
G.
,
Ge
,
L.
, et al
,
2021
, “
A Soft Neuroprosthetic Hand Providing Simultaneous Myoelectric Control and Tactile Feedback
,”
Nat. Biomed. Eng.
,
7
, p.
12
.
7.
Gladman
,
A. S.
,
Matsumoto
,
E. A.
,
Nuzzo
,
R. G.
,
Mahadevan
,
L.
, and
Lewis
,
J. A.
,
2016
, “
Biomimetic 4D Printing
,”
Nat. Mater.
,
15
(
4
), pp.
413
418
.
8.
Bendsøe
,
M. P.
, and
Kikuchi
,
N.
,
1988
, “
Generating Optimal Topologies in Structural Design Using a Homogenization Method
,”
Comput. Methods Appl. Mech. Eng.
,
71
(
2
), pp.
197
224
.
9.
Zhu
,
B.
,
Zhang
,
X.
,
Zhang
,
H.
,
Liang
,
J.
,
Zang
,
H.
,
Li
,
H.
, and
Wang
,
R.
,
2020
, “
Design of Compliant Mechanisms Using Continuum Topology Optimization: A Review
,”
Mech. Mach. Theory
,
143
, p.
103622
.
10.
James
,
K. A.
,
Kennedy
,
G. J.
, and
Martins
,
J. R.
,
2014
, “
Concurrent Aerostructural Topology Optimization of a Wing Box
,”
Comput. Struct.
,
134
, pp.
1
17
.
11.
Carstensen
,
J. V.
,
Lotfi
,
R.
,
Guest
,
J. K.
,
Chen
,
W.
, and
Schroers
,
J.
,
2015
, “
Topology Optimization of Cellular Materials With Maximized Energy Absorption
,”
Proceedings of the ASME 2015 International Design Engineering Technical Conferences
,
Boston, MA
,
Aug. 2–5
, pp.
1
10
.
12.
Maute
,
K.
,
Tkachuk
,
A.
,
Wu
,
J.
,
Qi
,
H. J.
,
Ding
,
Z.
, and
Dunn
,
M. L.
,
2015
, “
Level Set Topology Optimization of Printed Active Composites
,”
ASME J. Mech. Des.
,
137
(
11
), p.
111402
.
13.
Bhattacharyya
,
A.
, and
James
,
K.
,
2021
, “
Topology Optimization of Shape-Memory Polymer Structures With Programmable Morphology
,”
Struct. Multidiscipl. Optim.
,
63
(
1
), pp.
1863
1887
.
14.
Sossou
,
G.
,
Demoly
,
F.
,
Belkebir
,
H.
,
Qi
,
H. J.
,
Gomes
,
S.
, and
Montavon
,
G.
,
2019
, “
Design for 4D Printing: Modeling and Computation of Smart Materials Distributions
,”
Mater. Des.
,
181
, p.
108074
.
15.
Sun
,
X.
,
Yue
,
L.
,
Yu
,
L.
,
Shao
,
H.
,
Peng
,
X.
,
Zhou
,
K.
,
Demoly
,
F.
,
Zhao
,
R.
, and
Qi
,
H. J.
,
2021
, “
Machine Learning-Evolutionary Algorithm Enabled Design for 4D-Printed Active Composite Structures
,”
Adv. Funct. Mater.
,
32
(
10
), p.
2109805
.
16.
Kang
,
Z.
, and
James
,
K.
,
2019
, “
Multimaterial Topology Optimization for Elastic and Thermal Response
,”
Int. J. Numer. Methods Eng.
,
117
(
10
), pp.
1019
1037
.
17.
Haque
,
M. A.
,
Lavrik
,
N. V.
,
Hensley
,
D.
,
Briggs
,
D. P.
, and
McFarlane
,
N.
,
2021
, “
Carbonized Polymer for Joule Heating Processing Towards Biosensor Development
,”
Proceedings of the 2021 43rd Annual International Conference of the IEEE Engineering in Medicine & Biology Society (EMBC)
,
Mexico
,
Nov. 1–5
, pp.
7578
7581
.
18.
Paul
,
R. P.
,
1981
,
Robot Manipulators: Mathematics, Programming, and Control: The Computer Control of Robot Manipulators
,
MIT Press
,
Cambridge, MA
.
19.
Baghani
,
M.
,
Naghdabadi
,
R.
,
Arghavani
,
J.
, and
Sohrabpour
,
S.
,
2012
, “
A Thermodynamically-Consistent 3D Constitutive Model for Shape Memory Polymers
,”
Int. J. Plast.
,
35
, pp.
13
30
.
20.
Svanberg
,
K.
,
1987
, “
The Method of Moving Asymptotes—A New Method for Structural Optimization
,”
Int. J. Numer. Methods Eng.
,
24
(
2
), pp.
359
373
.
21.
Bendsøe
,
M. P.
, and
Sigmund
,
O.
,
1999
, “
Material Interpolation Schemes in Topology Optimization
,”
Arch. Appl. Mech.
,
69
(
9–10
), pp.
635
654
.
22.
Hu
,
M.
,
Ge
,
X.
,
Chen
,
X.
,
Mao
,
W.
,
Qian
,
X.
, and
Yuan
,
W.-E.
,
2020
, “
Micro/Nanorobot: A Promising Targeted Drug Delivery System
,”
Pharmaceutics
,
12
(
7
), p.
665
.
23.
Lipson
,
H.
, and
Kurman
,
M.
,
2013
,
Fabricated: The New World of 3D Printing
,
Wiley
,
Indianapolis, IN
.