Abstract

Avoiding stress concentrations is essential to achieve robust parts since failure tends to originate at such concentrations. With recent advances in multimaterial additive manufacturing, it is possible to alter the stress (or strain) distribution by adjusting the material properties in selected locations. Here, we investigate the use of grayscale digital light processing (g-DLP) 3D printing to create modulus gradients around areas of high stress. These gradients prevent failure by redistributing high stresses (or strains) to the neighboring material. The improved material distributions are calculated using finite element analysis. The much-enhanced properties are demonstrated experimentally for thin plates with circular, triangular, and elliptical holes. This work suggests that multimaterial additive manufacturing techniques like g-DLP printing provide a unique opportunity to create tougher engineering materials and parts.

1 Introduction

Avoiding failure through material design is an essential focus of the materials community for structural and functional applications. Several approaches have been proposed to achieve more robust materials, such as molecular design by incorporating multiple interactions [1,2], multiphase materials [35] to dissipate energy, and architectural materials using multimaterial fabrication [68]. Similarly, functionally graded materials (FGMs) with mechanical gradients have been used to satisfy various applications in aerospace, biomedicine, and defense [913]. FGMs with composition or stiffness variation in one system can also be found abundantly in nature. The heterogeneous structure and properties enable many functional features and enhanced mechanical properties [11,14]. For example, bone-to-tendon connections comprise a soft tendon and stiff bone joined by a strong interfacial bond, enabling significant energy absorption capabilities [15].

Using a mechanical property gradient is an effective method to tune the stress distribution and minimize failure in load­bearing structures. Stress concentrations occur in regions of a part where the geometry has caused stresses to build up during loading. This is crucial to the strength of a part since failure tends to initiate in the regions with the highest stress concentration factors (SCFs). The traditional method to reduce stress concentration is to adjust the geometry. A sharp corner, for example, can be replaced with a curved corner (or a fillet) to reduce the SCF. However, this process alters the part’s final shape, which may be undesirable in certain instances. An early approach to reduce the SCF around a hole was to introduce additional, smaller holes around the original one, but this method still alters the geometry and introduces additional, albeit lower, stress concentrations around the new holes [16]. Several other approaches have been developed to combat stress concentrations around defects without altering the geometry. In some cases, fibers with different orientations are used to provide different degrees of reinforcement at key points around the hole [1719]. Shah et al. used strategically placed piezoelectric patches to reduce stresses in certain regions [20]. Another approach has employed polymers with tunable stiffness to create radially symmetric gradients around circular holes to lower SCFs [2123]. Some research has already been conducted evaluating the effectiveness of grayscale printing to redistribute the stresses around defects by creating regions of locally varying Young’s modulus. Leben et al. [24] attempted to use an optimization approach in which material points throughout the plate were allowed to take a different modulus to reduce the strain energy gradient throughout the part. Additionally, Gu et al. [25] demonstrated the potential of using tailored stiffness distributions to improve resistance to crack propagation. While many of these approaches were successful, there is still much room for improvement. Specifically, fiber and piezo-patch reinforcement methods are tedious to fabricate and have a small design space. Methods such as radially symmetric stiffness gradients also have a limited design space and cannot be easily applied to geometries other than circular holes. Optimization or machine learning approaches have more versatility and much larger design spaces, but calculating a viable design could take a long time [24].

Over the past decade, additive manufacturing, or 3D printing (3DP), has enjoyed a steady rise in research and industrial realms with broad applications in metamaterials [2629], soft robotics [3032], high-performance materials [3336], and biomedical devices [3739]. The advent of 3DP technologies offers an attractive pathway toward the fabrication of FGMs by enabling more heterogeneous architectures and more complex mechanical gradients than conventional manufacturing processes [40,41]. This is enabled by the ability of 3DP to spatially control the local microstructure and the chemical composition in a layer-by-layer fashion. Currently, there are many successful attempts at printing graded materials with stiffness variation. The most common method is PolyJet printing, which uses rows of inkjet nozzles of the size of a few tens of microns to deposit inks of tailored compositions [42]. However, stringent ink requirements and the high cost of PolyJet printers limit their wide applications. Multinozzle direct ink writing (DIW) printing has also been widely used to fabricate FGMs with tunable mechanical property gradients [21,43,44]. Kokkinis et al. [21] reported a DIW-based multimaterial 3D­printing method to fabricate FGMs with modulus spanning three orders of magnitude and investigated bioinspired gradient designs for defect-tolerant structures. However, due to restrictions on the printing resolution of DIW, manipulating the mechanical properties in a continuous and high-resolution manner is very difficult. More recently, Kuang et al. [45] developed grayscale digital processing (g-DLP) printing using a two-stage curing ink that utilizes a thermal postcuring process to achieve widely tunable mechanical property gradients in a high-resolution manner. DLP printing offers several advantages over other printing methods. First, the process is generally faster as it cures each layer in one light exposure, which can be further sped up by using the continuous liquid interface production technique [46]. Second, DLP printers can use a variety of photocurable resins to achieve a wide range of material properties and dimension control [47,48]. Furthermore, through g-DLP printing, it is possible to produce parts containing a spectrum of material properties from a single resin with the potential to modulate the distribution of stresses generated within.

The present work proposes a generalized method to design modulus distributions that can significantly increase a part’s toughness compared to the equivalent uniform design. The proposed method locally adjusts the modulus of a part based on the intensity of the strain in the same location. The strain field is computed using finite element analysis (FEA). The entire design process is self-contained within the commercial FEA software package abaqus. With slight changes to key parameters, the method can be easily applied to different geometries or loading conditions. Coupling this robust design process with g-DLP printing allows the unique modulus distributions to be efficiently fabricated. We present material distribution designs with increased toughness and improved failure properties without altering the geometry of a part. This method allows for the production of complex components that possess features traditionally considered defects (e.g., holes and sharp corners) for applications in which the shape of the part cannot be compromised or adjusted.

2 Materials and Methods

2.1 Polymer Formulation.

The material used in this work is similar to the one described by Kuang et al. [45]. The ink contains 79.49 wt% bisphenol A ethoxylate diacrylate (avg. Mn 468) (BPAEDA), 19.87 wt% glycidyl methacrylate, 0.60 wt% bis(2,4,6-trimethylbenzoyl)-phenylphosphineoxide (PI 819) photoinitiator, and 0.04 wt% Sudan I dye. The aforementioned chemicals are purchased from Sigma Aldrich Inc. (St. Louis, MO, USA).

2.2 Grayscale 3D Printing.

The continuously graded samples are fabricated using g-DLP 3D printing. In the g-DLP printing process, the light patterns projected into the photocurable resin vat can vary spatially in intensity throughout the sample, unlike the uniform intensity patterns used in traditional DLP printing. Modulating the light intensity allows for the local control of the polymer network formation, leading to different material properties [49,50]. After a thermal postcuring stage, a single resin can achieve over three orders of magnitude in Young’s modulus [45].

Figure 1(a) shows a schematic of the bottom-up g-DLP printer used. It consists of a DLP PRO4500 projector (Wintech Digital System Technology Corp., Carlsbad, CA, USA) and an MTS50-Z8 (Thorlabs Inc, Newton, NJ, USA) linear motion stage. During printing, an image is projected toward the liquid resin, which is between a transparent window and the build platform or the previously printed layer. After a layer is cured, the linear stage shifts up to detach the current layer from the transparent window and allow the resin to flow into the void, to become the next layer. An example of a grayscale image is shown in Fig. 1(b). In the image, the brighter regions correspond to areas that will receive higher light intensities, whereas the darker regions receive less. The images are displayed one at a time, as seen in Fig. 1(c).

Fig. 1
Printing process via g-DLP printing: (a) a schematic of a bottom-up g-DLP printer, (b) an example of a grayscale image that can be printed, (c) layer order, (d) plot of light intensity versus modulus of printed parts, and (e) modulus of the printed part based on a grayscale image.
Fig. 1
Printing process via g-DLP printing: (a) a schematic of a bottom-up g-DLP printer, (b) an example of a grayscale image that can be printed, (c) layer order, (d) plot of light intensity versus modulus of printed parts, and (e) modulus of the printed part based on a grayscale image.
Close modal

2.3 Computational Methods.

This work uses the abaqus/standard fea package (Dassault Systèmes, Waltham, MA, USA) to calculate the numerical model predictions. To build a model with a spatially variable modulus, we create a user material model, which is based on the Neo-Hookean hyperelastic material with a locally varying modification factor multiplied to the shear modulus at each integration point. For each design, the model has displacement boundary conditions, and a three-step simulation is conducted. In the first step, the model is loaded (by displacement control), and the strain field is calculated. The modulus modification factor is 1.0 in the first loading step. The model is then unloaded to the initial state in step two. In step three, the load (displacement) is reapplied with modulus modification factors applied uniquely to each integration point, and the result is compared to step one.

2.4 Printed Sample Testing.

Both materials for characterization and samples from defect tolerance design are tested on a 10 kN electromechanical load cell (MTS, Eden Prairie, MN). Testing is performed at a strain rate of 1% per second. The models tested in abaqus are printed with solid regions on either end of the region of interest for mounting in the testing unit. The additional regions have a curved profile, flaring out away from the part. This shape is inspired by the standard dog bone tensile test design. The FEA models are altered to account for the additional material. Strain data are collected using a video extensometer setup during tensile tests. The video is processed through Adobe After Effects (Adobe Inc., San Jose, CA, USA), and the strain data are plotted with the stress data obtained from the tensile test by matching datapoints at corresponding time values.

3 Results

3.1 Mechanical Properties of g-DLP Printed Materials.

The mechanical properties of a photocurable resin change as a function of the degree of cure. This can be accomplished nonuniformly within a single part through g-DLP printing. For this article, the mechanical modulus is of notable interest. Figure 1(d) shows the modulus versus the light intensity. The light intensity refers to an integer value between 0 and 255, specifying the brightness of the projector pixel. A brightness of 255 corresponds to an in-plane intensity of 3.10 mW/cm2. The data for Fig. 1(d) are collected from testing dog bone samples, which are all printed at uniform intensities. All samples are exposed to UV light for the same duration of 16 s for each layer except the first which is cured for 35 s to improve adhesion to the build plate. As shown in Fig. 1(d), the resin experiences a tenfold change in modulus from 11 to 110 MPa over the 131–245 intensity range. This range of moduli is used for the remainder of the article. Figure 1(e) shows a modulus map of the part from the image in Fig. 1(b).

The stress–strain curves can be found in Appendix  A. In this work, strain is considered for the modulus modification process over stress because we assume constant failure strain independent of modulus. Failure cannot be inferred from the stress profile within a part because failure stress is a function of modulus, and by design, the elastic modulus will vary postmodification. It is therefore necessary to have a modulus-independent parameter. We support this assumption of constant failure strain with the data shown in Supplemental Fig. 1(b) available in the Supplemental Materials on the ASME Digital Collection., in which the materials printed from 150 to 255 fail at approximately the same strain independent of stiffness.

3.2 Modification Method.

To modify the modulus distribution, a baseline loading of the desired part is simulated first using displacement control to obtain the necessary modification factors at each integration point of the FEA model. For generality, the stiffness of the region of interest is taken as 1.0 MPa, considered the overall low value analogous to the lowest modulus applicable to the true material. The distribution of the maximum principal strain (ε1) at the end of the loading step is recorded. The strain values are mapped to the new modulus, such that regions of high-strain map to the upper limit of the elastic modulus range, while low-strain regions remain relatively unchanged (at the initial low stiffness). This remapping is performed using a modified sigmoid function to isolate key parameters,
f(Δε1,i)=1+A11+exp(B(Δε1,iC))
(1)
where A is the normalized maximum value of the property range (computed as the largest modulus in the property range divided by the smallest), B modulates the steepness of the curve, C affects the offset of the output, and Δε1,i represents the normalized strain at the ith integration point. The normalization of Δε1,i considers the overall minimum and maximum strain values in the part, computed as follows:
Δε1,i=ε1,iε1,minε1,maxε1,min
(2)
where ε1,max and ε1,min are the maximum and minimum values of the max principal strain across the part, respectively.

The outputs f(Δε1,i) of the mapping function are stored in a user-defined variable and then applied as the new modulus in the third step of the simulation, which deforms the sample in the same manner as the first step for comparison.

3.3 Modification Results.

We first present a simple case of a thin plate with a circular hole in horizontal tension. The test area is a 15 mm × 15 mm with a hole of 2 mm radius placed at the center. This geometry is ideal for the initial trials, given its simplicity and geometric symmetries. For this problem, one edge is fixed, and the opposite edge is displaced by 0.15 mm.

Figure 2 shows the process and performance of the model. The simulation and printed results are presented using values of A = 10, B = 15, and C = 0.45 in Eq. (1), and the corresponding plot of the sigmoid function is presented in Fig. 2(a). Figure 2(b) shows the distribution of ε1 in the unmodified model (only the square region of interest is shown for clarity). As expected, the largest strains occur above and below the hole. In contrast, the areas to the left and right of the hole show relatively small strains. The modulus is updated and reloaded using the strain pattern of Fig. 2(b). The corresponding strain distribution and modulus pattern are shown in Figs. 2(c) and 2(d), respectively. The maximum strain in Fig. 2(c) decreases by 53% in comparison to the unmodified sample in Fig. 2(b), and the overall distribution of strain is much closer to a uniform spread. All of this is promising for the prospect of increasing toughness.

Fig. 2
Design process of the modification algorithm: (a) plot of Eq. (1) using the specified B and C parameters for the tested sample, (b) ε1 distribution of the initial unmodified sample postloading depicted without the dog bone wings in view, (c) ε1 distribution of the modified postloading, (d) modulus pattern of the modified part depicted with wings, (e) pattern with uniform modulus with equivalent stiffness as the modified part, (f) grayscale DLP image of the modified part, (g) samples of untested print (top), tested modified print (middle), and tested uniform print (bottom), and (h) averaged stress–strain curves for the modified and uniform printed samples.
Fig. 2
Design process of the modification algorithm: (a) plot of Eq. (1) using the specified B and C parameters for the tested sample, (b) ε1 distribution of the initial unmodified sample postloading depicted without the dog bone wings in view, (c) ε1 distribution of the modified postloading, (d) modulus pattern of the modified part depicted with wings, (e) pattern with uniform modulus with equivalent stiffness as the modified part, (f) grayscale DLP image of the modified part, (g) samples of untested print (top), tested modified print (middle), and tested uniform print (bottom), and (h) averaged stress–strain curves for the modified and uniform printed samples.
Close modal

A control sample is printed to determine if the modifications improve the performance of the part. To obtain the most direct comparison, the control sample is assigned a uniform modulus across the 15 mm × 15 mm area such that the two samples have an equivalent effective stiffness. This equivalent stiffness is calculated from the displacement from edge to edge of the whole 15 mm × 15 mm area of the modified sample (which needs to be calculated because as the part becomes stiffer, the wings deform slightly), and the reaction force generated from the displacement load. As such, each uniform control sample will have a unique stiffness dependent on the parameters of the modified sample. Figure 2(e) shows the corresponding uniform sample with equal stiffness to the modified sample of Fig. 2(d). The grayscale image for the modified sample is provided in Fig. 2(f). The required grayscale values are determined from the property versus intensity relationship shown in Fig. 1(d).

Figure 2(g) shows printed g-DLP samples. The top two samples are those with the modified modulus patterns (before and after being tested). The bottom print is a tested sample of the uniform pattern. Despite the distinct difference between the strain patterns in Figs. 2(b) and 2(c) (the pattern in Fig. 2(b) can be assumed to be the same for any sample of the uniform modulus), the modified sample still fails at the same region as the uniform sample, demonstrating the dominance of the stress concentration formed by the hole. The stress–strain curves in Fig. 2(h) are obtained by averaging the stress–strain curves of three specimens in both groups. The unprocessed stress–strain curves for the tensile tests can be found in Appendix  B. The method of averaging stress is shown in Appendix  C. As shown in Fig. 2(h), despite failing in similar regions, the overall performance of the modified and uniform samples is not equivalent. The modified samples show an almost 90% increase in toughness (computed as the area under the average stress–strain curve) and a 35% increase in failure strain as compared to the uniform ones. The improvement in toughness is likely the result of the softer regions of the modified part. The soft regions to the left and right of the hole, shown in Fig. 2(d), deform before the stiffer regions above and below the hole, causing the highest strains to develop away from the hole first. Since there are lower strains near the defect, larger loads are required to cause the part to fail. The overlap between the stress–strain curves indicates excellent agreement between their overall stiffnesses.

3.4 Optimization.

The results in Fig. 2 show that the proposed methodology has the potential to improve the failure strain of a part without compromising its stiffness. However, the arbitrary selection of parameters for Eq. (1) suggests the need for further optimization. To do this, we perform a property sweep across a range of possible remapping parameters to identify higher- and lower-performing designs.

To quantitatively rank the improvement of parameters across the property sweep, a metric analyzing the final strain distribution of a modified part is needed. It was found that using the entire field to compute, a metric is generally unnecessary. The proposed metric instead only considers the maximum value of ε1, which we assume is where a part is most likely to fail. Also, decreasing/increasing ε1 at the primary defect will increase/decrease the likelihood of failure. The values of the metric are summarized in the following equation:
M=ε1,maxuε1,maxmε1,maxu
(3)

The superscript “m” denotes the modified model, and “u” is the corresponding uniform pattern. The index “max” refers to the integration point containing the maximum value of ε1 found in the uniform model. Both ε1,maxu and ε1,maxm are considering the same integration point/location even though this index may not correspond to the integration point of maximum ε1 in the modified model. Every term in Eq. (3) is a function of A, B, and C and the construction of the model (e.g., geometry, boundary conditions, and loading). With this metric, a positive value corresponds to a decrease in ε1 from uniform to modified and a negative value corresponds to an increase.

The following property sweep only considers an A value of 10 because it corresponds to the fixed printable property range shown in Fig. 1(d). Ten different B values in the range of −10 to 60 and ten values of C in the range of 0.1 to 0.55 are considered culminating in 100 simulations. Figure 3 details the results.

Fig. 3
Optimization study of three different samples, (a) property sweep of the proposed metric, (b) plot of Eq. (1) for the B and C parameters of the “inferior,” “neutral,” and “improved” samples tested, (c) modulus pattern of the “inferior” sample, (d) modulus pattern of the “neutral” sample, (e) modulus pattern of the “improved” sample, (f) ε1 distribution of the “inferior” modified sample postloading, (g) ε1 distribution of the “neutral” modified sample postloading, (h) ε1 distribution of the “improved” modified sample postloading, (i–k) samples of the modified (top) and uniform (bottom) prints, (i) printed “inferior” designs, (j) printed “neutral” designs, (k) printed “improved” designs, (l) video frame during failure of an “improved” modified sample, and (m) averaged stress–strain curves for the modified and uniform printed samples for all three designs.
Fig. 3
Optimization study of three different samples, (a) property sweep of the proposed metric, (b) plot of Eq. (1) for the B and C parameters of the “inferior,” “neutral,” and “improved” samples tested, (c) modulus pattern of the “inferior” sample, (d) modulus pattern of the “neutral” sample, (e) modulus pattern of the “improved” sample, (f) ε1 distribution of the “inferior” modified sample postloading, (g) ε1 distribution of the “neutral” modified sample postloading, (h) ε1 distribution of the “improved” modified sample postloading, (i–k) samples of the modified (top) and uniform (bottom) prints, (i) printed “inferior” designs, (j) printed “neutral” designs, (k) printed “improved” designs, (l) video frame during failure of an “improved” modified sample, and (m) averaged stress–strain curves for the modified and uniform printed samples for all three designs.
Close modal

The distribution of the described metric is plotted in Fig. 3(a). The regions (B < 0) correspond to samples with increased ε1 at the defect, the regions (B ≥ 0, C ≤ ~0.3) correspond to samples of little to no overall change, and the green regions (B ≥ 0, C ≥ ~0.3) correspond to samples of moderate to high improvement (decreased ε1 at the defect). The maximum value of the plot is above 0.8, which indicates an 80% decrease in maximum strain from uniform to modified. Three samples are chosen from different regions of the property sweep to confirm the utility of our design metric. These samples are referred to as “inferior” (B = −10, C = −0.35), “neutral” (B = 25, C = 0.10), and “improved” (B = 60, C = 0.35) corresponding to the regions indicated by the circle, diamond, and star, respectively. The corresponding sigmoid functions are detailed in Fig. 3(b). Notice that choosing a negative B value flips the sigmoid curve about the midpoint of the y-axis (in this case, the line drawn at f(Δε1,i)=5.5), resulting in stiffening areas of low strain while not altering areas of high strain.

Figures 3(c)3(e) show the modulus patterns of the inferior, neutral, and improved samples postmodification. As expected, the distributions of Figs. 3(c) and 3(e) appear almost inverted. Figures 3(f)3(h) show the ε1 distributions of the inferior, neutral, and improved samples postmodification and reloading. Observe that the inferior sample shown in Fig. 3(f) appears like an exaggerated version of the strain distribution of a uniform sample, as shown in Fig. 2(b), which suggests a very early failure. The visual improvements in Fig. 3(g) appear minimal compared to a uniform sample, most noticeably the high-strain regions above and below the holes, suggesting a low change in toughness. Lastly, the improved sample possesses a radically different strain distribution shown in Fig. 3(h), the hole showing little to no strain concentrations with most of the strain concentrated at the far edges of the sample.

Figures 3(i)3(k) show printed samples of the modified (top images) and uniform (bottom images) inferior, neutral, and improved samples after testing. For Figs. 3(i) and 3(j), the samples fail at the hole, similar to what is shown in Fig. 2(g). However, the failure in the improved sample of Fig. 3(k) is very different, failing away from the hole. Figure 3(l) is a frame from a video taken while testing one of the improved samples. One can see the crack beginning to propagate through the sample. The crack does not form at the interface between the modified region and the wings (which are visible in the frame), instead emerging in one of the four regions of high strain near the far edges of the part shown in Fig. 3(h). The shape of the deformed hole is also highly deformed with the top and bottom sides remaining almost undeformed, while the left and right sides experience substantial strain.

Figure 3(m) shows the averaged stress–strain extensometer curves for the inferior, neutral, and improved modified and uniform samples. The shapes of the modified and their respective uniform samples have an excellent agreement. The failure of all the uniform samples is at a strain of 0.15 ± 0.01. The trends produced by the change in failure strain of each modified sample compared to its respective uniform sample match with the expectation described in the metric (the improved samples behave better than the inferior samples). However, the precise change in failure strain is not accurately predicted by the metric as the metric stated an 81% decrease for the inferior sample, but the failure strain only decreases by about 36%. Similarly, the metric predicts an 85% increase for the improved sample, but the failure strain increases by 66%. The neutral samples are in excellent agreement between metric and change in failure strain. The metric can predict general trends toward improvements but not accurately predict exact changes. However, this is expected since the metric is based on a singular number over a wide distribution. Additionally, the FEA model lacks failure modeling, so the metric extrapolates part failure from an elastic load.

A special two-hole model is constructed to demonstrate the application of this process with deterministic failure. The two-hole model is constructed of two 7.5 mm wide by 15 mm tall regions with holes in their centers, attached side by side to form a 15 mm × 15 mm plate. We intend to show that a plate with multiple defects can be tuned to have a controlled failure defect without compromising the stiffness of the whole structure. In this model, the right region is modified with the desired pattern from given A, B, and C parameters, while the left region is given an equivalent uniform stiffness to match the right. The process for generating the modulus pattern is the same as described in Sec. 3.3, except during the initial loading, unloading, and reloading with the modified right region, the left region is given the same stiffness as the wings to minimize its impact on the right region. One additional step is added to the simulation in which the left region is updated to its new uniform value and the whole part is reloaded.

Figure 4 shows two different patterns for the two-hole model. An “inferior” and “improved” (analogous to what is described in Fig. 3) sample are chosen to show different options for controlled failure. The property sweep shown in Fig. 3(a) is used to guide the choice of patterns. Figures 4(a) and 4(b) show the modulus patterns for the two samples chosen. Figures 4(c) and 4(d) show one print for either pattern. Three samples of each are printed and tested. The right halves of the printed samples possess the modified region, while the left halves have the appropriately modified stiffnesses. The inferior samples always failed on the modified (right) side, and the improved samples always failed on the uniform (left) side. These further demonstrate the consistent failure properties of the modified regions.

Fig. 4
Two-hole sample modification study: (a) modulus pattern of the “inferior” design, (b) modulus pattern of the “improved” design, (c) one printed and tested sample of the “inferior” design. The hole on the left is in the uniform region, and the hole on the right is in the modified region and (d) one printed and tested sample of the “improved” design. The hole on the left is in the uniform region, and the hole on the right is in the modified region.
Fig. 4
Two-hole sample modification study: (a) modulus pattern of the “inferior” design, (b) modulus pattern of the “improved” design, (c) one printed and tested sample of the “inferior” design. The hole on the left is in the uniform region, and the hole on the right is in the modified region and (d) one printed and tested sample of the “improved” design. The hole on the left is in the uniform region, and the hole on the right is in the modified region.
Close modal

3.5 Samples With Equilateral Triangle Holes.

Thus far, the only geometry considered has been a simple plate with a circular hole cut out of the center. It is important to consider other geometries to demonstrate the generality of the proposed modification algorithm. However, in the interest of maintaining simplicity and obtaining results that can be reasonably compared with the previous examples, the new geometries will continue the style of a plate in tension with a hole through the center. The hole shape becomes the new design parameter in the following models. Here, we first consider the sample with equilateral triangle holes, as shown in Fig. 5.

Fig. 5
Modification and optimization study of a plate with a triangular hole: (a) property sweep of the triangular hole model, (b) ε1 distribution of the initial unmodified sample postloading, (c) modulus pattern of “improved” design postmodification, (d) modulus pattern of “neutral” design postmodification, (e) averaged stress–strain curves for the modified and uniform printed samples for the two designs, (f) ε1 distribution of the “improved” design postmodification and reloading, (g) printed modified “improved” design postloading, (h) Two-hole variant of triangular hole model using the same design parameters as the “improved” design, and (i) Two printed samples of the two-hole “improved” design.
Fig. 5
Modification and optimization study of a plate with a triangular hole: (a) property sweep of the triangular hole model, (b) ε1 distribution of the initial unmodified sample postloading, (c) modulus pattern of “improved” design postmodification, (d) modulus pattern of “neutral” design postmodification, (e) averaged stress–strain curves for the modified and uniform printed samples for the two designs, (f) ε1 distribution of the “improved” design postmodification and reloading, (g) printed modified “improved” design postloading, (h) Two-hole variant of triangular hole model using the same design parameters as the “improved” design, and (i) Two printed samples of the two-hole “improved” design.
Close modal

A triangularly shaped hole is selected because it is the simplest shape with corners and straight edges and can be oriented so that it is only symmetric about one axis. In this test, the triangle is oriented so that it “points” in the direction of loading (i.e., symmetric about the x-axis only). However, the corners of an equilateral triangle would create enormous stress concentrations and skew the modification results. To correct this, the corners are all filleted with a fillet radius equal to 5% of the final height of the triangle. The overall height of the triangular hole is 4 mm, with fillets on each corner with 0.2 mm radii.

A new property sweep is needed since the geometry of the model has now changed. For this new property sweep, six values of B in the range of −10 to 25 and eight values of C in the range of 0.1 to 0.45 were used. The A parameter remains at 10. In preliminary simulations, the resolution and breadth of the modulus pattern are more sensitive to the B parameter since the high stress/strain concentrations at the corners of the hole produced large differences between the minimum and maximum strain values. At higher B values, only the regions directly adjacent to the corners receive any degree of modification, while the rest of the sample is not altered, resulting in very soft parts. The range of B values is decreased since it becomes impractical to print and test such soft samples. The new property sweep is shown in Fig. 5(a).

Figure 5(b) shows the distribution of ε1 for the initial unmodified sample. As mentioned earlier regarding the sensitivity to the B parameter, the highest strain regions occur in a small zone adjacent to the corners. This small region would translate to only a few pixels of the grayscale image emitting any high amount of UV light from the projector.

For comparison and relating to the property sweep, two patterns are selected, one that falls in the “improved” (green, B > 0, C > 0.2) region and the other in the “neutral” (blue, B ≥ 0, C ≤ 0.2) region of Fig. 5(a). The modulus patterns for the improved (B = 35, C = 0.25) and neutral (B = 5, C = 0.10) patterns are shown in Figs. 5(c) and 5(d), respectively. The samples are only picked from the improved and neutral regions because there is a reasonable expectation of a noticeable difference in their failure strain, but they are not intentionally designed to fail prematurely like the inferior samples. The improved and neutral regions also possess B values that are strictly greater than zero, making their corresponding sigmoid functions more comparable.

Figure 5(e) shows the average stress–strain curves of the improved and neutral designs alongside their equivalent uniform designs. There is a close agreement between the shape of the modified and uniform curves, similar to the results shown in Fig. 3(m). The improved and neutral samples perform as well as expected based on the trends in the property sweep. The neutral modified samples fail at nearly the same point as the corresponding uniform samples on average despite the property sweep predicting an approximate 20% increase in failure strain from the uniform to the modified. In contrast, the improved samples experience an approximate 81% increase in failure strain, while the value obtained from the metric predicted a 76% increase.

To visualize the impact of the improved modulus pattern, Fig. 5(f) shows the distribution of ε1 post reloading. In this strain distribution, the highest strain values appear far away from the hole, decreasing the strain values at the corners of the hole to under 25% of the initial case seen in Fig. 5(b). However, despite the toughness increase of the improved sample, it still fails at the corners, as shown in Fig. 5(g).

As a further demonstration of the superior performance of the improved sample, a two-hole model is constructed with the same triangular hole using the same sigmoid function. The modulus pattern of the modified sample is presented in Fig. 5(h), and two instances of the printed version are shown in Fig. 5(i). As expected in the tests, the two-hole model consistently fails at the uniform hole.

3.6 Samples With Rotated Ellipse Holes.

The second of the geometries tested is a rotated ellipse (Fig. 6). An ellipse is selected because it is a smooth shape without any corners while having clear stress concentration points. The hole was rotated to make it asymmetric about both the x and y-axes. The ellipse has a three-to-one aspect ratio between its major and minor axes and is rotated 30 deg clockwise from the y-axis. The major axis of the ellipse has a length of 7.5 mm, and the minor axis measures 2.5 mm.

Fig. 6
Modification and optimization study of a plate with a rotated elliptical hole: (a) property sweep of the rotated elliptical hole model, (b) ε1 distribution of the initial unmodified sample postloading, (c) modulus pattern of “improved” design postmodification, (d) modulus pattern of “neutral” design postmodification, (e) averaged stress–strain curves for the modified and uniform printed samples for the two designs, (f) ε1 distribution of the “improved” design postmodification and reloading, (g) printed modified “improved” design postloading, (h) two-hole variant of the rotated elliptical hole model using the same design parameters as the “improved” design, and (i) two printed samples of the two-hole “improved” design.
Fig. 6
Modification and optimization study of a plate with a rotated elliptical hole: (a) property sweep of the rotated elliptical hole model, (b) ε1 distribution of the initial unmodified sample postloading, (c) modulus pattern of “improved” design postmodification, (d) modulus pattern of “neutral” design postmodification, (e) averaged stress–strain curves for the modified and uniform printed samples for the two designs, (f) ε1 distribution of the “improved” design postmodification and reloading, (g) printed modified “improved” design postloading, (h) two-hole variant of the rotated elliptical hole model using the same design parameters as the “improved” design, and (i) two printed samples of the two-hole “improved” design.
Close modal

The same A, B, and C parameters as the triangular hole model are explored for the property sweep of this model for the same reasons presented above. The property sweep is presented in Fig. 6(a). The elliptical hole model seems less sensitive to the negative B values than the triangular and circular hole models.

Figure 6(b) shows the ε1 distribution for the initial unmodified sample postloading, which experiences similar issues of large strains in very small regions like the triangular hole sample. The modulus pattern of the improved (B = 35, C = 0.30) and neutral (B = 10, C = 0.20) samples selected are shown in Figs. 6(c) and 6(d), respectively. Comparing the improved and neutral modulus patterns, the advantages of conservatively reinforcing select regions of the part are clear, mainly around defects, as it allows other regions the flexibility to absorb the applied load. Figure 6(e) shows the average stress–strain curves of the improved and neutral models chosen alongside their equivalent uniform models. The improved model experiences an 85% increase in failure strain, while the predicted change is 76%. Additionally, the neutral model experiences a 38% increase with a predicted change of 34%.

Figure 6(f) shows the ε1 distribution for the improved pattern. The regions containing the highest values of ε1 shift to the side of the stress concentration points of the ellipse. Figure 6(g) shows an example of one of the modified, improved samples posttesting. This image suggests the sample fails slightly away from the peak of the ellipse, as predicted from the strain distribution of Fig. 6(f). However, this occurrence is inconsistent, as modified samples occasionally fail at the peaks of the ellipse.

A two-hole variant of the elliptical hole model is presented as a final demonstration. Figure 6(h) shows the modulus pattern for this model, using the same A, B, and C parameters used in the improved sample. Figure 6(i) shows the two printed samples failing on the uniform side, demonstrating improved agreement with the predicted outcome.

4 Conclusion

A new approach is presented for producing functionally graded materials via g-DLP printing with improved mechanical properties (most notably increased toughness). Using the maximum principal strain distributions produced by FEA, the proposed method considers stiffening areas of a part under high strain while keeping low-strain areas soft. Postmodification, the improved designs demonstrate a reduction in strain at their primary defects and an increase in strain in regions previously subject to small strains. These new regions experiencing high strain exhibit low elastic stiffness and possess little to no stress-concentrating geometry, making the whole part less liable to failure. The modification procedure is predicated on the assumption that the failure strain of g-DLP printing parts is independent of the degree of UV curing observed in experiments. A property sweep over the tunable parameters is performed to investigate the impact of changing such parameters on the improvement of the part and determine if an optimal modification exists. A metric is proposed that condensed the overall improvement of a model down to a single number. The proposed metric is found to be able to predict trends in failure strain. Despite the choice of g-DLP printing, the proposed modification algorithm should be generalizable to other manufacturing methods capable of locally controlling the modulus over a continuous range of values. However, the performance metric may need to be adjusted based on the failure properties of the materials used. Lastly, three different hole geometries are simulated and printed to demonstrate the generalizability of the proposed methodology. Altogether, the proposed modification procedure shows the ability to produce parts with remarkable improvements in mechanical toughness and large decreases in stress concentrations without alterations to the part’s geometry or stiffness.

Acknowledgment

The support of an AFOSR grant (A9550-20-1-0306; Dr. B.-L. “Les” Lee, Program Manager) and an NSF grant (1931977) are gratefully acknowledged. Gift funds from HP, Inc and Northrop Grumman Corporation are also greatly appreciated. This article has been authored by an employee of National Technology & Engineering Solutions of Sandia, LLC under Contract No. DE-NA0003525 with the U.S. Department of Energy (DOE). The employee owns all right, title, and interest in and to the article and is solely responsible for its contents. The United States Government retains and the publisher, by accepting the article for publication, acknowledges that the United States Government retains a nonexclusive, paid-up, irrevocable, world-wide license to publish or reproduce the published form of this article or allow others to do so, for United States Government purposes. The DOE will provide public access to these results of federally sponsored research in accordance with the DOE Public Access Plan.2

Footnote

Conflict of Interest

There are no conflicts of interest.

Data Availability Statement

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.

Appendix A: Tensile Test Results of Dog Bone Samples

A key assumption of the proposed modification process is that the cured resin fails at the same strain independent of the light dose. To justify this assumption, dog bone samples are printed at varying light intensities. The failure strength and stiffness are computed through a standard tensile test. The stiffnesses obtained are then used to construct the curve in Fig. 1(c), the curve terminates at 245 because that produces a clean tenfold stiffness change from a value of 11 MPa to 110 MPa. Five light intensities are tested in the range of 130–255. A cure time of 35 s is used for the first layer and 16 s for each subsequent layer.

As shown in Supplemental Fig. 1(a), the samples tested from 150 to 255 all exhibit similar failure strains consistent with the assumption made, but below 150, the failure strains decrease as the light intensity decreases. However, this diminishing failure strain has a little impact on the results compiled in the modification. None of the samples with a uniform modulus distribution are printed with an intensity around 130 or 140, as such their failure strains are within the constant failure strain region (between 150 and 255). Most of the modified parts do contain areas exposed to light in the 130, 140 range; however, by design, such regions experience low-strain loads due to the geometric conditions of the model and are unlikely to contribute to the failure. This reasoning is further justified by the results shown in Fig. 3(m), in which all the uniform samples fail at similar strains, and the modified “neutral” part failed at the same strain as the corresponding uniform sample.

The exception is the inferior sample shown in Fig. 3, this sample contains regions of high strain (above and below the hole) that receives a very low light dose. However, it is unlikely this contributes significantly to the failure strain since the strain at the hole computed in the FEA (Fig. 3(f)) is over twice that computed in the neutral sample FEA (Fig. 3(g)) implying it will fail sooner. Additionally, the improved samples shown in Fig. 3 do not fail at the hole, instead failing in a region receiving a very low light dose, and this suggests that if the primary assumption about uniform failure strain were true, and these samples may perform even better.

Appendix B: Tensile Test of Uniform and Modified Samples Corresponding to Fig. 2 

Unaveraged stress–strain data from the tensile tests for the samples tested in Fig. 2 are presented. Due to their higher strength, the modified samples tend to fail abruptly, as is shown by the sudden drop-off in stress in Supplemental Fig. 2(a). With this sudden failure, the samples will initially fail either on the top or bottom (top-to-bottom aligned with the y-axis, orthogonal to the direction of loading) of the hole, which will allow them to recover before fully failing on the other side. This recovery is visible on test 4 of Supplemental Fig. 2(a), and in this trial, the sample fails at a strain of 0.2 and then again at a strain of 0.22. However, for consistency, only the first failure event is considered. As such for test 4, all data after a strain of 0.2 are truncated in regard to the analysis.

In contrast, the uniform samples would often fail more gradually, due to their relative softness at the primary defect. This is most apparent in test 3 of Supplemental Fig. 2(b), the stress appears to almost plateau at a strain of 0.12, before fully failing at a strain of 0.14. This occurs because the part gradually rips apart at the strain of 0.12, with cracks on either side of the hole slowly propagating outward simultaneously. To properly analyze this failure, all data after the start of failure (a strain of 0.12 for test 3) are ignored.

Appendix C: Averaging Method of Tensile Test Curves

The average curves shown in Figs. 2(h), 3(m), 5(e), and 6(e) are computed from fitted exponential curves based on the stress–strain data points collected from the tensile tests. The data points are fitted to curves because it allows for simple interpolation, permitting a direct comparison of stresses from different tests at the same strains since the data obtained from the tensile test and video extensometer are not at consistent strain increments across experiments. The fitted curves are computed using the fit function in matlab’s curve fitting toolbox only using the data collected before the initial failure of the respective test. All fitted curves are prescribed in the following form:
σ¯i=aiexp(biε)+ciexp(diε)
where σ¯i is the fitted stress in MPa of the ith dataset in a sample group (e.g., all the curves in Supplemental Fig. 2(a) form one group and the stress and strain data for “Modified, Test 1” is one dataset), ε is the engineering strain, and ai, bi, ci, and di are the fitted constants computed in fit. This two-exponential function is selected because it presents consistent agreement across all datasets collected from the tensile test. The domain used for ε across a sample group is between zero and the arithmetic average of all the failure strains in that group. Another advantage of the fitted curve is that it can produce a continuation of the dataset beyond the failure strain obtained in the experiment, allowing a dataset to contribute to the average even if the strain value is beyond the failure strain for that given experiment. The final averaged curve, σ¯ave, computed over the n datasets of a group is
σ¯ave(ε)=1ni=1nσ¯i(ε)

References

1.
Zhang
,
Y. S.
, and
Khademhosseini
,
A.
,
2017
, “
Advances in Engineering Hydrogels
,”
Science
,
356
(
6337
), p.
eaaf3627
.
2.
Sun
,
T. L.
,
Kurokawa
,
T.
,
Kuroda
,
S.
,
Ihsan
,
A. B.
,
Akasaki
,
T.
,
Sato
,
K.
,
Haque
,
M. A.
,
Nakajima
,
T.
, and
Gong
,
J. P.
,
2013
, “
Physical Hydrogels Composed of Polyampholytes Demonstrate High Toughness and Viscoelasticity
,”
Nat. Mater.
,
12
(
10
), pp.
932
937
.
3.
Liff
,
S. M.
,
Kumar
,
N.
, and
McKinley
,
G. H.
,
2007
, “
High-Performance Elastomeric Nanocomposites via Solvent-Exchange Processing
,”
Nat. Mater.
,
6
(
1
), pp.
76
83
.
4.
Pei
,
A.
,
Malho
,
J.-M.
,
Ruokolainen
,
J.
,
Zhou
,
Q.
, and
Berglund
,
L. A.
,
2011
, “
Strong Nanocomposite Reinforcement Effects in Polyurethane Elastomer With Low Volume Fraction of Cellulose Nanocrystals
,”
Macromolecules
,
44
(
11
), pp.
4422
4427
.
5.
Kuang
,
X.
,
Chen
,
K.
,
Dunn
,
C. K.
,
Wu
,
J.
,
Li
,
V. C. F.
, and
Qi
,
H. J.
,
2018
, “
3D Printing of Highly Stretchable, Shape-Memory and Self-Healing Elastomer Toward Novel 4D Printing
,”
ACS Appl. Mater. Interfaces
,
10
(
8
), pp.
7381
7388
.
6.
Wang
,
Z.
,
Xiang
,
C.
,
Yao
,
X.
,
Le Floch
,
P.
,
Mendez
,
J.
, and
Suo
,
Z.
,
2019
, “
Stretchable Materials of High Toughness and Low Hysteresis
,”
Proc. Natl. Acad. Sci.
,
116
(
13
), pp.
5967
5972
.
7.
Libonati
,
F.
, and
Buehler
,
M. J.
,
2017
, “
Advanced Structural Materials by Bioinspiration
,”
Adv. Eng. Mater.
,
19
(
5
), p.
1600787
.
8.
Dimas
,
L. S.
,
Bratzel
,
G. H.
,
Eylon
,
I.
, and
Buehler
,
M. J.
,
2013
, “
Tough Composites Inspired by Mineralized Natural Materials: Computation, 3D Printing, and Testing
,”
Adv. Funct. Mater.
,
23
(
36
), pp.
4629
4638
.
9.
Sola
,
A.
,
Bellucci
,
D.
, and
Cannillo
,
V.
,
2016
, “
Functionally Graded Materials for Orthopedic Applications—An Update on Design and Manufacturing
,”
Biotechnol. Adv.
,
34
(
5
), pp.
504
531
.
10.
Gupta
,
A.
, and
Talha
,
M.
,
2015
, “
Recent Development in Modeling and Analysis of Functionally Graded Materials and Structures
,”
Prog. Aerosp. Sci.
,
79
, pp.
1
14
.
11.
Chen
,
P.-Y.
,
McKittrick
,
J.
, and
Meyers
,
M. A.
,
2012
, “
Biological Materials: Functional Adaptations and Bioinspired Designs
,”
Prog. Mater. Sci.
,
57
(
8
), pp.
1492
1704
.
12.
Liu
,
Z.
,
Meyers
,
M. A.
,
Zhang
,
Z.
, and
Ritchie
,
R. O.
,
2017
, “
Functional Gradients and Heterogeneities in Biological Materials: Design Principles, Functions, and Bioinspired Applications
,”
Prog. Mater. Sci.
,
88
, pp.
467
498
.
13.
Roach
,
D. J.
,
Hamel
,
C. M.
,
Wu
,
J.
,
Kuang
,
X.
,
Dunn
,
M. L.
, and
Qi
,
H. J.
,
2017
, “
4-D Printing: Potential Applications of 3-D Printed Active Composite Materials
,”
J. Homeland Defense Secur. Inf. Anal. Center
,
4
(
4
), pp.
20
27
.
14.
Wegst
,
U. G. K.
,
Bai
,
H.
,
Saiz
,
E.
,
Tomsia
,
A. P.
, and
Ritchie
,
R. O.
,
2015
, “
Bioinspired Structural Materials
,”
Nat. Mater.
,
14
(
1
), pp.
23
36
.
15.
Rossetti
,
L.
,
Kuntz
,
L. A.
,
Kunold
,
E.
,
Schock
,
J.
,
Müller
,
K. W.
,
Grabmayr
,
H.
,
Stolberg-Stolberg
,
J.
, et al
,
2017
, “
The Microstructure and Micromechanics of the Tendon–Bone Insertion
,”
Nat. Mater.
,
16
(
6
), pp.
664
670
.
16.
Jindal
,
U. C.
,
1983
, “
Reduction of Stress Concentration Around a Hole in a Uniaxially Loaded Plate
,”
J. Strain Anal. Eng. Des.
,
18
(
2
), pp.
135
141
.
17.
Muc
,
A.
, and
Ulatowska
,
A.
,
2012
, “
Local Fibre Reinforcement of Holes in Composite Multilayered Plates
,”
Compos. Struct.
,
94
(
4
), pp.
1413
1419
.
18.
Wang
,
G.
,
2018
, “
Homogenized and Localized Stress Reconfigurations of Solid or Hollow Fiber Reinforced Materials in a Multi-Scale Framework
,”
Compos. Struct.
,
184
, pp.
1099
1110
.
19.
Boddeti
,
N.
,
Tang
,
Y.
,
Maute
,
K.
,
Rosen
,
D. W.
, and
Dunn
,
M. L.
,
2020
, “
Optimal Design and Manufacture of Variable Stiffness Laminated Continuous Fiber Reinforced Composites
,”
Sci. Rep.
,
10
(
1
), p.
16507
.
20.
Shah
,
D. K.
,
Joshi
,
S. P.
, and
Chan
,
W. S.
,
1994
, “
Stress Concentration Reduction in a Plate With a Hole Using Piezoceramic Layers
,”
Smart Mater. Struct.
,
3
(
3
), pp.
302
308
.
21.
Kokkinis
,
D.
,
Bouville
,
F.
, and
Studart
,
A. R.
,
2018
, “
3D Printing of Materials With Tunable Failure via Bioinspired Mechanical Gradients
,”
Adv. Mater.
,
30
(
19
), p.
1705808
.
22.
Yang
,
Q.
,
Gao
,
C.-F.
, and
Chen
,
W.
,
2010
, “
Stress Analysis of a Functional Graded Material Plate With a Circular Hole
,”
Arch. Appl. Mech.
,
80
(
8
), pp.
895
907
.
23.
Sburlati
,
R.
,
2013
, “
Stress Concentration Factor Due to a Functionally Graded Ring Around a Hole in an Isotropic Plate
,”
Int. J. Solids Struct.
,
50
(
22
), pp.
3649
3658
.
24.
Leben
,
L. M.
,
Schwartz
,
J. J.
,
Boydston
,
A. J.
,
D’Mello
,
R. J.
, and
Waas
,
A. M.
,
2018
, “
Optimized Heterogeneous Plates With Holes Using 3D Printing via Vat Photo-Polymerization
,”
Addit. Manuf.
,
24
, pp.
210
216
.
25.
Gu
,
G. X.
,
Chen
,
C.-T.
,
Richmond
,
D. J.
, and
Buehler
,
M. J.
,
2018
, “
Bioinspired Hierarchical Composite Design Using Machine Learning: Simulation, Additive Manufacturing, and Experiment
,”
Mater. Horiz.
,
5
(
5
), pp.
939
945
.
26.
Zheng
,
X.
,
Lee
,
H.
,
Weisgraber
,
T. H.
,
Shusteff
,
M.
,
DeOtte
,
J.
,
Duoss
,
E. B.
,
Kuntz
,
J. D.
, et al
,
2014
, “
Ultralight, Ultrastiff Mechanical Metamaterials
,”
Science
,
344
(
6190
), pp.
1373
1377
.
27.
Zhang
,
H.
,
Guo
,
X.
,
Wu
,
J.
,
Fang
,
D.
, and
Zhang
,
Y.
,
2018
, “
Soft Mechanical Metamaterials With Unusual Swelling Behavior and Tunable Stress-Strain Curves
,”
Sci. Adv.
,
4
(
6
), p.
eaar8535
.
28.
Wu
,
S.
,
Eichenberger
,
J.
,
Dai
,
J.
,
Chang
,
Y.
,
Ghalichechian
,
N.
, and
Zhao
,
R. R.
,
2022
, “
Magnetically Actuated Reconfigurable Metamaterials as Conformal Electromagnetic Filters
,”
Adv. Intell. Syst.
,
4
(
9
), p.
2200106
.
29.
Ma
,
C. P.
,
Wu
,
S.
,
Ze
,
Q.
,
Kuang
,
X.
,
Zhang
,
R.
,
Qi
,
H. J.
, and
Zhao
,
R.
,
2021
, “
Magnetic Multimaterial Printing for Multimodal Shape Transformation With Tunable Properties and Shiftable Mechanical Behaviors
,”
ACS Appl. Mater. Interfaces
,
13
(
11
), pp.
12639
12648
.
30.
Kim
,
Y.
,
Yuk
,
H.
,
Zhao
,
R.
,
Chester
,
S. A.
, and
Zhao
,
X.
,
2018
, “
Printing Ferromagnetic Domains for Untethered Fast-Transforming Soft Materials
,”
Nature
,
558
(
7709
), pp.
274
279
.
31.
Rus
,
D.
, and
Tolley
,
M. T.
,
2015
, “
Design, Fabrication and Control of Soft Robots
,”
Nature
,
521
(
7553
), pp.
467
475
.
32.
Roach
,
D. J.
,
Hamel
,
C. M.
,
Dunn
,
C. K.
,
Johnson
,
M. V.
,
Kuang
,
X.
, and
Qi
,
H. J.
,
2019
, “
The m4 3D Printer: A Multi-Material Multi-Method Additive Manufacturing Platform for Future 3D Printed Structures
,”
Addit. Manuf.
,
29
, p.
100819
.
33.
Lee
,
C.-U.
,
Vandenbrande
,
J.
,
Goetz
,
A. E.
,
Ganter
,
M. A.
,
Storti
,
D. W.
, and
Boydston
,
A. J.
,
2019
, “
Room Temperature Extrusion 3D Printing of Polyether Ether Ketone Using a Stimuli-Responsive Binder
,”
Addit. Manuf.
,
28
, pp.
430
438
.
34.
Hegde
,
M.
,
Meenakshisundaram
,
V.
,
Chartrain
,
N.
,
Sekhar
,
S.
,
Tafti
,
D.
,
Williams
,
C. B.
, and
Long
,
T. E.
,
2017
, “
3D Printing All-Aromatic Polyimides Using Mask-Projection Stereolithography: Processing the Nonprocessable
,”
Adv. Mater.
,
29
(
31
), p.
1701240
.
35.
Kuang
,
X.
,
Zhao
,
Z.
,
Chen
,
K.
,
Fang
,
D.
,
Kang
,
G.
, and
Qi
,
H. J.
,
2018
, “
High-Speed 3D Printing of High-Performance Thermosetting Polymers via Two-Stage Curing
,”
Macromol. Rapid Commun.
,
39
(
7
), p.
1700809
.
36.
Matsuzaki
,
R.
,
Ueda
,
M.
,
Namiki
,
M.
,
Jeong
,
T.-K.
,
Asahara
,
H.
,
Horiguchi
,
K.
,
Nakamura
,
T.
,
Todoroki
,
A.
, and
Hirano
,
Y.
,
2016
, “
Three-Dimensional Printing of Continuous-Fiber Composites by In-Nozzle Impregnation
,”
Sci. Rep.
,
6
(
1
), p.
23058
.
37.
Colosi
,
C.
,
Shin
,
S. R.
,
Manoharan
,
V.
,
Massa
,
S.
,
Costantini
,
M.
,
Barbetta
,
A.
,
Dokmeci
,
M. R.
,
Dentini
,
M.
, and
Khademhosseini
,
A.
,
2016
, “
Microfluidic Bioprinting of Heterogeneous 3D Tissue Constructs Using Low-Viscosity Bioink
,”
Adv. Mater.
,
28
(
4
), pp.
677
684
.
38.
Qiu
,
K.
,
Zhao
,
Z.
,
Haghiashtiani
,
G.
,
Guo
,
S.-Z.
,
He
,
M.
,
Su
,
R.
,
Zhu
,
Z.
, et al
,
2018
, “
3D Printed Organ Models With Physical Properties of Tissue and Integrated Sensors
,”
Adv. Mater. Technol.
,
3
(
3
), p.
1700235
.
39.
Derby
,
B.
,
2012
, “
Printing and Prototyping of Tissues and Scaffolds
,”
Science
,
338
(
6109
), pp.
921
926
.
40.
Studart
,
A. R.
,
2016
, “
Additive Manufacturing of Biologically-Inspired Materials
,”
Chem. Soc. Rev.
,
45
(
2
), pp.
359
376
.
41.
Porter
,
M. M.
,
Ravikumar
,
N.
,
Barthelat
,
F.
, and
Martini
,
R.
,
2017
, “
3D-Printing and Mechanics of Bio-Inspired Articulated and Multi-Material Structures
,”
J. Mech. Behav. Biomed. Mater.
,
73
(
Supplement C
), pp.
114
126
.
42.
Cazón
,
A.
,
Morer
,
P.
, and
Matey
,
L.
,
2014
, “
PolyJet Technology for Product Prototyping: Tensile Strength and Surface Roughness Properties
,”
Proc. Inst. Mech. Eng. B
,
228
(
12
), pp.
1664
1675
.
43.
Chen
,
K.
,
Zhang
,
L.
,
Kuang
,
X.
,
Li
,
V.
,
Lei
,
M.
,
Kang
,
G.
,
Wang
,
Z. L.
, and
Qi
,
H. J.
,
2019
, “
Dynamic Photomask-Assisted Direct Ink Writing Multimaterial for Multilevel Triboelectric Nanogenerator
,”
Adv. Funct. Mater.
,
29
(
33
), p.
1903568
.
44.
Hansen
,
C. J.
,
Saksena
,
R.
,
Kolesky
,
D. B.
,
Vericella
,
J. J.
,
Kranz
,
S. J.
,
Muldowney
,
G. P.
,
Christensen
,
K. T.
, and
Lewis
,
J. A.
,
2013
, “
High-Throughput Printing via Microvascular Multinozzle Arrays
,”
Adv. Mater.
,
25
(
1
), pp.
96
102
.
45.
Kuang
,
X.
,
Wu
,
J.
,
Chen
,
K.
,
Zhao
,
Z.
,
Ding
,
Z.
,
Hu
,
F.
,
Fang
,
D.
, and
Qi
,
H. J.
,
2019
, “
Grayscale Digital Light Processing 3D Printing for Highly Functionally Graded Materials
,”
Sci. Adv.
,
5
(
5
), p.
eaav5790
.
46.
Tumbleston
,
J. R.
,
Shirvanyants
,
D.
,
Ermoshkin
,
N.
,
Janusziewicz
,
R.
,
Johnson
,
A. R.
,
Kelly
,
D.
,
Chen
,
K.
, et al
,
2015
, “
Continuous Liquid Interface Production of 3D Objects
,”
Science
,
347
(
6228
), pp.
1349
1352
.
47.
Montgomery
,
S. M.
,
Demoly
,
F.
,
Zhou
,
K.
, and
Qi
,
H. J.
,
2023
, “
Pixel-Level Grayscale Manipulation to Improve Accuracy in Digital Light Processing 3D Printing
,”
Adv. Funct. Mater.
, p.
2213252
.
48.
Tanaka
,
M.
,
Montgomery
,
S. M.
,
Yue
,
L.
,
Wei
,
Y.
,
Song
,
Y.
,
Nomura
,
T.
, and
Qi
,
H. J.
,
2023
, “
Turing Pattern-Based Design and Fabrication of Inflatable Shape-Morphing Structures
,”
Sci. Adv.
,
9
(
6
), p.
eade438
.
49.
Zhao
,
Z.
,
Wu
,
J.
,
Mu
,
X.
,
Chen
,
H.
,
Qi
,
H. J.
, and
Fang
,
D.
,
2017
, “
Origami by Frontal Photopolymerization
,”
Sci. Adv.
,
3
(
4
), p.
e1602326
.
50.
Zhao
,
Z.
,
Wu
,
J.
,
Mu
,
X.
,
Chen
,
H.
,
Qi
,
H. J.
, and
Fang
,
D.
,
2017
, “
Desolvation Induced Origami of Photocurable Polymers by Digit Light Processing
,”
Macromol. Rapid Commun.
,
38
(
13
), p.
1600625
.

Supplementary data