Multimaterial 4D printing with a tunable bending model

Shape-memory polymer (SMP)-based functional structures may now be produced more efficiently via four-dimensional (4D) printing, benefiting from the recent advances in multi-material three-dimensional printing technologies. Composite material design using 4D printing has opened new possibilities for customizing the shape memory property of smart polymers. This work studies a design strategy to harness desirable morphing by 4D printing multimaterial composites with a focus on the detailed finite element (FE) procedure, experimental results, and soft robotic application. Composites with bilayer laminates consisting of a SMP and a flexible elastomer are constructed with variable thickness ratios to control the self-bending of the composite. FE simulations are used to understand the underlying processes of composite materials and to generate accurate predictions for the experimental results, which reduces cost and development time. The application of 4D printing and multi-material composite programming is demonstrated with a soft robotic gripper for manipulating fragile objects.


Introduction
With the advent of four-dimensional (4D) printing, it has become much simpler to create shape memory polymer (SMP)-based functional structures. Among the many types of active materials, SMPs are among the most common because of their ability to achieve a programmed shape memory * Authors to whom any correspondence should be addressed.
Original content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. effect in response to external stimuli [1]. Actuators [2], flexible electronics [3], soft robotics [4], biomedical devices [5], aerospace [6], and others [7][8][9][10][11] have all been shown to benefit from SMPs because of their stiffness variable characteristic and efficient processing methods.
In practice, heat responsive SMPs have proven to be the most feasible method for these applications using 4D printing. Wu et al [12] developed 4D-printed multi-shape active composites by alternating the position of fibers in the top and bottom layer, instigating on the specific advantages of additive manufacturing compared to traditional composite laminates assembling. Yuan et al [13] studied the tunable shape memory behavior of 4D-printed bilayer laminates of heat responsive SMPs and elastomers in different layouts concerning their shape fixity to obtain bending motions. Teoh et al [14] used two basic materials provided by Stratasys (USA), namely VeroWhitePlus and TangoBlackPlus, for 4D printing of selffolding flowers in an ambient environment.
One-dimensional composite rods were designed and manufactured to yield desired three-dimensional (3D) structures in a reproducible approach by Ding et al, in which the authors demonstrate the ability of 4D printing to vary the elastomer layer height fraction and the cross-section twist angle along the rod axis to self-assemble into complex 3D shapes [15]. Bodaghi et al [16] employed self-morphing 4D printing as actuation elements to design complex structures that show planar to 3D shape-shifting. The effects of printing speed on the self-morphing characteristics have been investigated in their work. Alsheby et al used various bioinspired geometrical shapes and patterns in 4D printing of heat responsive soft actuators [17]. Shapes of different sizes and combinations of different patterns were investigated to achieve different curvature control of 4D-printed structures. The fundamental requirements of 4D printing modeling, simulation, and hands on practice on heat responsive SMPs were demonstrated by Zolfagharian et al [18]. Pandini et al [19] studied the sequential self-deployment of auxetic structures utilizing the 4D printing technology. Gries et al [20] discovered applications for heat responsive 4D printing in smart textiles. More research on remotely induced heat responsive 4D printing, such as magnetic, electric, and light stimuli, has recently been reported in order to broaden the applications of this emerging technology [21][22][23].
Composite material design using 4D printing has opened new possibilities for customizing the shape memory property. However, research on the potential of shape programming and reversibility in 4D printing is mostly qualitative at this time [24]. Typically, the morphed shapes after applying the stimulation were approximations based on the initial design, with no quantitative analysis. It is difficult to develop useful 4Dprinted devices due to the lack of predictability since the parameters of the transformed shape are left undetermined.
In this work, a design strategy to harness desirable morphing by 4D printing multi-material composites is investigated. Composites with bilayer laminates consisting of SMP and flexible elastomer are constructed in variable thickness ratios to control the bending of the composite. Finite element (FE) simulations in ABAQUS are used to study and anticipate experimental findings to understand the underlying processes of composite materials and to generate accurate prediction for the experimental data. In terms of FE model, the main novelty of this project is that the bending angle [18] of strips has been predicted in FE software. These strips have been made with two types of material (VeroBlack and Tan-goPlus), material properties are various, and it will change with temperature, further, the thickness of these materials in all of configuration is different. In this study, a comprehensive FE model has been presented that it can predict the bending angle of composite strips in various conditions. The application of 4D printing and multi-material composite programming is demonstrated with a soft robotic gripper for manipulating fragile objects. With its detailed and efficient modeling strategy and the practical demonstration, several areas of 4D printing, including bioprinting, soft robotics, vibration isolation, and self-morphing structures, might benefit from this method.
The remainder of this article is organized as follows. Section 2 describes the material properties, fabrication methods and the procedure for the FE simulations, which includes a detailed description of model geometry and properties, boundary conditions and mesh sizes in ABAQUS. The experimental and finite element model (FEM) results are discussed in section 3, which also includes a demonstration of the experimental testing and thermomechanical behavior for various layout designs with varying material compositions. Section 4 concludes this article and outlines the future work.

Materials and manufacturing
The design strategy proposed here is based on strain mismatches between the elastomer and the glassy polymer layers of a 3D-printed layered composite. Using a commercial Stratasys J750 3D printer and commercial materials Vero-Black (glassy polymer) and TangoPlus (elastomer), composite beams are 3D printed, which are then activated by heating in a uniform constant temperature.
The conceptual design is described here. A compressive stress/strain (pre-stress or pre-strain) is introduced into the elastomer, but not the glassy polymer, during the 3D printing process when the elastomer and glassy polymer are printed into a layered composite. At room temperature the elastomer is sufficiently constrained by the glassy polymer, which has a much higher Young's modulus (1 GPa versus 1 MPa) and a good bonding with the elastomer. Assuming the layers bond perfectly from a mechanical load distribution standpoint, the composite form will not change when the support is removed, and a high-fidelity and stable shape will be printed. However, heating the bilayer composite sample results in two specific phenomena, see figure 1(a). Initially, when the glassy polymer is heated over its glass transition temperature (T g ), its Young's modulus lowers dramatically, and the internal stress/strain built up in the off-center elastomer is released. This off-center compressive stress could be large enough to bend the bilayer and make the composite beam into an arch form. Second, the thermal strain mismatch becomes more prominent when the two materials are heated, since the elastomer has a greater coefficient of thermal expansion compared to the glassy polymer. The synergic effects of main drive coming from the stress/strain release in the programmed pre-strained/stressed elastomer, thermal expansions plus the dimensional factor [13] result in a self-bending feature for the bilayer composite beam. The shape recovery ratio of the materials, defined as the ratio between the recovered shape and the pre-programmed shape, is expected to be 100%. Once the composite beam is cooled down, well below T g of the glassy material, the glassy layer becomes hard and keep the arch form in the presence of compressive stress in the elastomer. The difference between the height of arch at the end of heating and (a) 4D-printed flat composite strip (the temporary shape) and then, when heated, transform into a curved shape (the permanent shape) that remains largely unchanged upon further cooling and heating cycles, the black and blue represents Vero and Tango, respectively. (b) Representative of thermomechanical and shape memory characteristics of the formative materials. The left column displays the temperature-dependent strain response for VeroBlack and TangoPlus throughout a complete shape memory cycle. The right column displays the dynamic mechanical analysis (DMA) for both materials. The results in (b) are reproduced from [13] under a Creative Commons Attribution (CC BY) license by Elsevier. Reproduced with permission from [13]. © 2020 Published by Elsevier Ltd. CC BY-NC-ND 4.0.
cooling steps can be associated to the shape fixity ratio. Later, it will be found that this ratio is large as high as 95%-98% that means the composite beam keeps 95%-98% of the maximum height at the end of cooling. Upon re-heating the composite beam, it returns to full 100% arch form. It implies that the beam will keep its arch form with just 2%-5% of changes during heating-cooling cycles.
Polymerized from isobornyl acrylate, acrylic monomer, urethane acrylate, epoxy acrylate, acrylic monomer, acrylic oligomer, and a photo initiator, VeroBlack exhibits the rigidity of plastic at room temperature. Polymerized from urethane acrylate oligomer, exo-1,7,7-trimethylbicyclo hept-2-yl acrylate, methacrylate oligomer, polyurethane resin, and photo initiator, TangoPlus has the flexibility of rubber. The printing process, which involves droplet deposition spreading into a continuous layer and photocuring, sets underneath a residual compressive stress/strain in the elastomer (see figure 1(a)). The compressive strain was determined to be 4.4% for the elastomer and 0% for the glassy polymer by means of photography, as described in [13]. Figure 1(b) displays VeroBlack's temperature-dependent strain response throughout the course of a complete shape memory cycle. The glass transition phenomenon, in which polymers undergo a reversible transition between a glassy state and a rubbery state, may be used to explain the shape memory behavior shown in figure 1(b) from the standpoint of polymer physics. As demonstrated in figure 1(b), the storage modulus of VeroBlack rises by nearly two orders of magnitude when the material is cooled from the rubbery state to the glassy one. By increasing the temperature from 20 to 80 degrees Celsius, the polymer chains are made very mobile, allowing the material to return to its undeformed structure on a laboratory time scale. Unlike VeroBlack, however, TangoPlus fully recovers from the strain it experiences at 80 • C when unloaded at 20 • C ( figure 1(b) bottom left). Since T g (1 • C) of TangoPlus is well below the programming temperature range (20 • C-80 • C), the material remains in its rubbery state even after cooling to room temperature, 20 • C, and the instantaneous shape recovery after unloading is presented by the high polymer chain mobility and entropic elasticity.
Based on this knowledge of the basic characteristics of the composite materials, five composites with varying thickness ratios but the same volume were created. The model is made by two materials and the thickness of Vero and Tango are different in each sample. The 4D-printed structure in this work can be investigated considering the following design parameters: length (L), width (W), and thickness (h 1 and h 2 ) of the shape memory printed polymer, as shown in figure 2 and table 1. The employed PolyJet multimaterial 3D printer (J750, Stratasys, Edina, MN, USA) can accurately inject multiple materials with an in-plane precision of 200 µm. The STL files were prepared and transferred into the PolyJet Studio software for the Stratasys J750 printer, where they were printed at a slice height of 27 µm and a resolution of 600 × 300 dots per inches in plane using the device's default high mix printing mode. Nozzles on the 3D printer were kept at 65 degrees Celsius throughout the production process. A layer printing time of 82 s was used for all samples. The geometrical information, shape memory characterization and the thermomechanical properties of the polymers such as Young's modulus, glass transition temperature, and the coefficient of thermal expansion, were measured according to [11] using rectangular strips samples with dimensions of 15 mm × 3 mm × 0.5 mm. The results are shown in tables 1-3.

Governing equations of the thermal-mechanical coupling
This section describes the detailed governing equations used for the composite 4D printing model in FE software. In this study, a new approach for the simulation of these composite strips is introduced. This method is based on changes in thermal expansion coefficient and degree of pre-strain. This model considers both parameters together (thermal expansion coefficient and pre-strain) and an effective thermal expansion coefficient is presented that it can predict the bending angle of actuators. The effective thermal expansion coefficient reflects the behavior of strips with changes in temperature. This model is not very complicated, and it has good accuracy. Furthermore, in this model, changes of materials properties with temperature have been considered, so it has a good accuracy in high temperature and after transition temperature (T g ). The CAD geometry for 3D-printed prototypes is imported to ABAQUS for simulation. In this case, boundary conditions and sources are predetermined, and a solid-state, physical heat transfer method is applied. The material library is consulted and updated in response to the criteria. The initial temperature of the samples is set at 298.15 K, which are then are heated to 343.15 K in an oven with constant temperature and the Vero half laying on the oven surface.
In the simulation of thermally sensitive 4D printing, mechanical stress is coupled with thermal load. Hooke's law for an isotropic material [25] in the context of thermo-elastic problems states the following: where ∆T, β, ε, D, σ are the gradient of the temperature profile, the thermal stress modulus, the strain, the matrix of elastic constants, thermal coefficient and the second term of equation points out to the Cauchy stress. a and β are determined as: where E Young's modulus, ν, Poisson's ratio, and α is the thermal expansion coefficient.
For a Cartesian 3D model, we may write the transient heat conduction equation as: where the thermal conductivities in the x, y, and z axes as a function of temperature are denoted by k x (T), k y (T) and k z (T). In addition, G is the heat generation per unit volume, T, temperature, t, time, ρ, density, C P , specific heat, and ε V , volume strain. To properly define the real-world problem, it is important to take into account both the existing boundaries and the initial condition. In addition, we are interested in the thermoelastic deformations that take place between 0 and t. The subsequent process involves establishing the matrix equations that fully characterize each element's characteristics. Left hand side (LHS [k]) matrices and load vectors ([f ]) are involved in these equations. After determining the LHS and load vector for each element, and defining [T] as global unknown vector, the global matrix for the whole area may be assembled as follows: In this equation, n represents the total number of nodes in an element, N i , represents the shape functions, and T i , represents the temperature at each node. The Galerkin technique is one of many possible approaches to solve this problem, and it is the one used here. An expression for the Galerkin technique used with equation (4) is: where Ω and Γ denote the domain of integration and the domain border, respectively. Furthermore, l, m, and n are the cosines of the corresponding outward surface normals. It is also assumed that the material is isotropic, leading to thermal conductivities k x = k y = k z . Applying the temperature relation in equation (6) into equation (7) results: where i and j stand for the x and y axes of nodes. The final equation (8) form is as follows: where [C ij ] =ˆρC P N i N j dΩ (10)

FE simulations
The FE software solution diagram shown in figure 3. In the following, the procedure for the simulation of the 4D-printed structures in ABAQUS is described.

Design the model geometry.
Here, we create a 3D model with a solid element within the part module. Figure 4 displays the tools needed to design geometric shapes, where the geometry is divided into two sections using the partition cell tool with the respective thickness for the Tango and Vero layers.

Set the model properties.
The core of the simulation is provided here. The properties of the Tango and Vero materials are updated in the property module, as listed in table 2 [13,26]. Then, each part is given the relevant characteristics. A change in temperature may cause a change in a property. Because transient heating is involved in the model, we need to specify the relationship between these characteristics and temperature using the Young's Modulus and the thermal expansion coefficient. Figure 5 displays the temperature dependence of the Young's Modulus for Vero and Tango [13]. For this stage, determining the coefficient of thermal expansion for each material is crucial. A positive thermal expansion coefficient is present in the material, but there is already some pre-strain present (negative expansioncontraction). The negative thermal expansion coefficient is reflected in the degree of pre-strain. This parameter is a combination of thermal expansion coefficient and pre-strain. These parameters are variable with temperature changes. In this numerical study, an effective thermal expansion coefficient is introduced that reflects the interaction of these parameters together. The effective thermal expansion coefficient is determined through a methodical approach that involved a systematic process to predict pre-strain values based on available data and previous knowledge. Table 3 shows effective thermal expansion coefficient for both materials as a function of temperature. The context of this module in ABAQUS is shown in figure 6.   [13,26].
Reference [13] 0.14 1050 2300 Vero Specific heat Cp J (Kg K) −1 6 1.4 1200 1300 To simulate this condition, a surface is considered that the samples are placed on it where the temperature of this surface is identical to oven's temperature. Figure 7 shows the environment of the assembly module.

Choosing a suitable solver.
Given that the mechanical load is a byproduct of the heating process in this investigation, a solution that takes into account both mechanical load and heat transfer is necessary. A coupled temperaturedisplacement scheme is the most effective method for solving this issue. This is a nonlinear calculation that considers the interaction of displacement and temperature simultaneously. For this purpose, an exact implementation of Newton's method involves a nonsymmetric Jacobian matrix as the coupled equations [27]. As can be seen in figure 8, a transient solution is chosen which accounts for temperature variability over time.

Define interactions.
As mentioned in the assembly section, the samples have been laid with the Vero side on the surface. The type of contact between the sample and the surface is normal behavior, which is also assumed frictionless. Furthermore, the surface-to-surface standard interaction has been chosen for this problem, as shown in figure 9.

Setting the loads and boundary conditions.
In this work, the main form of loading is temperature. The type of heat transfer is conduction by contact to the constant temperature surface (343.15 K). Also, the initial temperature for the samples is defined as 298.15 K. Moreover, the surface's boundaries are fixed (constrained) and two end edges of the samples are constrained. Figure 10 demonstrates the loads and boundary conditions for this problem.

Mesh study.
In FE simulations, mesh size is one of the most important parameters to improve convergence.
Most importantly, reliable results should be relatively meshindependent [28]. Hence, in this work, the bending angle of the 4D-printed structure is chosen as a parameter to evaluate the quality of the mesh. The bending angle is first determined by starting with a large seed size and gradually decreasing it until there is no longer any correlation between the two. This procedure is demonstrated in figure 11. As the only geometry variable parameter is thickness, Sample-1 has been chosen for the mesh sensitivity, because it has the least thickness among the other configurations and the size of its element is generalized for the other samples. The size of the seeds for both materials is the same.
The approximate size of the seed for the 4D-printed part in this simulation is 0.0012 m. There are 2135 nodes and 1440 elements for this part. The type of element is C3D8T. Additionally, the size of the seeds for the surface is 0.0025 m. There are 533 nodes and 480 elements for this part, and the element type is R3D4. The environment of the mesh module is shown in figure 12. An eight-node thermally connected brick with a linear temperature and displacement response is decided. The element is a hexagonal type.

Results and discussions
All the morphology studies using composite beams were conducted in an oven. Samples were placed in the oven preheated to the temperature 343.15 K. Between 20 and 40 s, they settled into their final shape. After 50 s, we captured the deformed configurations and used an image analysis software (ImageJ) to calculate the curvature of the transformed 4D-printed composites, ensuring the shape transformation was complete. The  maximum deformation happened at the end of the heating step. The shape recovery ratio of the materials, defined as the ratio between the recovered shape and the pre-programmed shape, was found 100% for all cases.
The predicted results from the FEM simulations are compared to the experimental results in figures 13-15, table 4 and  table 5. It is evident that the FEM predictions show good correlation with the experimentally observed shapes (see supplementary video 1). The main observation is the descending trend of the bending angle with larger Tango ratios in the composite beams, as seen in Samples 1-3, respectively. It is noteworthy that, even with the same ratio of composite materials, the bending angles decreased with higher total thickness, as seen in samples 4-6. Besides, table 4 reports the small error between the FEM and experimental results for the bending angles of the composite beams. Nevertheless, the FEM predictions are affected in samples 3 and 6 where the ration of Tango to Vero thickness stands more than double and the total thickness of the composite beam increases significantly, respectively.
Curvature morphologies with time for the 4D-printed composite beams are shown in figure 16, where the stated curvatures are the average of three series of measurements. For all results, the beams start to deform rapidly after being placed in the oven or with temperature rise. Subsequently, the sheet deforms along an upward convex direction and converges to a constant curvature. The time to achieve this shape changes differ for the beams as per the thickness ratio and total thickness of the composite materials. Samples 1-3 show higher bending rates in the initial stages, respectively. While samples 4-6 experience bending rates with descending speed, respectively, attributed to the greater total thickness. Also, all samples reach the maximum curvature before 120 s. In terms of the relatively large error for sample-3 and sample-6, there are different issues that could lead to this error. The manufacturing process and the rate of pre-strain depend on this process. The relatively large errors in sample-3 and sample-6 come from variables in the manufacturing process, and they can be improved by monitoring the conditions during the printing process. These findings indicate that as soon as the active layer temperature meets the glass transition point, the shrinking process initiates. The temperature gradient throughout the thickness of the composite beam causes the whole structure to curve and shrinking occurs at the equilibrium point between the passive and active layers. The adhesion interaction and the mechanical balance between active and passive parts may both have contributed to the minor disparities between the FEM and experimental results. All results presented so far have been related to the configuration of the composite beams at the end of the heating step. We also checked the configuration (height) of the beams after cooling to the room temperature and found that they keep 95%-98% of their height which is related to the shape fixity rate. 95%-98% for the shape fixity ratio is very satisfying and promising. It was also found that any further heating-cooling cycle leads the composite beam to switch a bit between full bending height and 95%-98% of it that means just 2%-5% height change during further heating-cooling cycles.
To demonstrate the application of multi-material 4Dprinted structures and the design procedure presented in this study, an ad hoc experiment in soft robotics was carried out. As shown in figure 17, two composite beams were 4D-printed as the soft grippers of a DOBOT Magician for handling delicate objects, such as lab glassware weighing 11 g. A blow dryer is used to generate the heat required for the bending of the 4D-printed gripper beams.

Conclusions
In this study, we showed how direct 4D printing could be used in the product design space for composite materials with self-bending features. As a proof of concept for controlling the bending, we built composites using bilayer laminates made of SMP and flexible elastomer with various thickness ratios. To better comprehend the mechanisms of action in composite materials and to provide an accurate forecast for the experimental data, FE model simulations were developed. The methodology, in essence, proposed a methodical design and production procedure based on the integrated assembly of different sections. Instead of a single configuration with known form and geometry in advance but little morphing capacity, one can then plan for two (or more) configurations' worth of geometry, materials, and manufacturing processes using the method presented in this paper. Many design opportunities are available due to the interplay of geometrical, material, and process characteristics. In this work, 2D flat beams were created and then turned into a 3D structure via 4D printing by using a purposefully designed material architecture. Although 2D beams has a basic geometry, the conceptual design and method supplied in this paper can be implemented to 4D print complicated structures in the future. The FEM developed in this model predicted the experimental results very well. As demonstrated in the results, the bending angle was decreased by increasing the overall thickness, and it changed between 40 and 70 degrees. Also, table 4 compared the results for experimental and FE. As seen in this table, the peak error is approximately 10%; this error is associated with the sample that has the highest thickness. Also, figure 16 shows that the samples will reach their final position between 20 and 40 s, and this behavior was very well estimated by the FEM.
The streamlined method of production of high-resolution complex 3D reprogrammable structures presented here holds great promise for a wide range of applications in fields as diverse as medicine, aerospace, and consumer goods. This also suggests a new paradigm in product design in which components are conceived to simultaneously embody multiple configurations during service. Inverse design or finding the starting 2D form and composition that will result in the desired 3D shape, is still an ongoing topic in 4D printing. Recently, machine learning and data-driven techniques have also emerged as an effective tool for resolving inverse design issues that may be augmented in this study for developing highly sophisticated, morphing, and architecturally designed materials.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).