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).
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.
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.
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].
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].
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 at the designated node, N (marked with a red dot in Fig. 5, at the end of the thermomechanical cycle (), while constraining the optimizer to ensure that the total volume of active material 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
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.
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.
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.
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.
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.
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.
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 , which produced a twist angle of per element as shown in Fig. 8(a). For the bending elements, eight different samples were tested with strains ranging between and . These results were then interpolated to obtain the expected bending angle for a strain of , which was the applied strain of the self-tying knot during thermomechanical programming. This interpolation produced a bending angle of for each individual bending element. These results ( and ) 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 and 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.
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.
For each element in the sequence, offset distances along each axis (δx, δy, and δz) are calculated for the displaced kinematic element.
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
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) */
/* Evaluate change in position vector */
Re = RzRyRx /* Evaluate the total rotation based on the current element using Eq. (10) */
R ← RRe /* Update the total rotation matrix R */
/* 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 a − g correspond to the letters of the seven kinematic element options shown in Fig. 4.
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.
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.
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).
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
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