Multi‐Material 3D Printing of Functionally Graded Hierarchical Soft–Hard Composites

Hard biological tissues (e.g., nacre and bone) have evolved for millions of years, enabling them to overcome the conflict between different mechanical properties. The key to their success lies in the combination of limited material ingredients (i.e., hard and soft constituents) and mechanistic ingredients (e.g., functional gradients and building block hierarchical organization). However, the contribution of each material and mechanistic ingredient is still unknown, hindering the development of efficient synthetic composites. Quantitative and systematic studies of hard–soft composites are required to unravel every factor's role in properties outcome. Herein, a voxel‐by‐voxel multi‐material 3D printing technique is used to design and additively manufacture different groups of hard–soft composites. Several combinations of gradients, multilevel hierarchies, and brick‐and‐mortar arrangements are created. Single‐edge notched fracture specimens are mechanically tested and computationally simulated using extended finite element method (XFEM). It is found that functional gradients alone are not sufficient to improve fracture properties. However, up to twice the fracture energy of the hard face is observed when combining functional gradients with hierarchical designs, significantly increasing composite properties. Microscopic analysis, digital image correlation, and strain distributions predicted with XFEM are used to discuss the mechanisms responsible for the observed behaviors.


Introduction
Hard biological tissues, of which bone and nacre are prime examples, are staggeringly efficient structural materials. Through millions of years of evolutionary refinement, nature has found ways to overcome the conflict between different mechanical properties. [1,2] The structural performances of nearly all conventional engineering materials pale in comparison with that high level of multifaceted mechanical efficiency. A particularly tall order that engineering materials usually fall short of is combining high levels of stiffness with toughness and strength. Engineers, therefore, have to choose one property over another, often resulting in toughness being prioritized over strength for the reasons of safety, [3,4] thereby rendering engineering structures highly inefficient. [1] It is unclear how exactly hard biological tissues achieve such remarkable degrees of structural performance. Intensive research, particularly during the last decade, has successfully identified the ingredients of this success story. On the one hand, the material ingredients are usually a relatively soft organic phase, which is based on natural polymers such as proteins and polysaccharides, and an inorganic hard phase, which is made from ceramics such as calcium salts or silica. [2,[5][6][7] On the other hand, the main design strategies or the mechanistic ingredients are the application of functional gradients, [8,9] and intricate placement of hard-soft structures (e.g., in a brick-and-mortar fashion) as well as the use of hierarchical structures. [10][11][12][13] However, identifying such material and mechanistic ingredients has not enabled us to achieve a similar level of efficiency in our man-made materials, as it is unclear how these mechanistic and material ingredients interact with each other at different spatial and time scales, how important the contribution of each mechanism is, and how one mechanism modulates the effects of the others. In a way, the major mysteries in the design of hard tissues remain largely unsolved.
One approach used by many researchers to unravel the mystery of hard tissues is detailed histological studies aimed at better understanding the nature and role of different material ingredients [14][15][16][17] as well as detailed mechanical characterization DOI: 10.1002/adem.201901142 Hard biological tissues (e.g., nacre and bone) have evolved for millions of years, enabling them to overcome the conflict between different mechanical properties. The key to their success lies in the combination of limited material ingredients (i.e., hard and soft constituents) and mechanistic ingredients (e.g., functional gradients and building block hierarchical organization). However, the contribution of each material and mechanistic ingredient is still unknown, hindering the development of efficient synthetic composites. Quantitative and systematic studies of hard-soft composites are required to unravel every factor's role in properties outcome. Herein, a voxel-by-voxel multi-material 3D printing technique is used to design and additively manufacture different groups of hard-soft composites. Several combinations of gradients, multilevel hierarchies, and brickand-mortar arrangements are created. Single-edge notched fracture specimens are mechanically tested and computationally simulated using extended finite element method (XFEM). It is found that functional gradients alone are not sufficient to improve fracture properties. However, up to twice the fracture energy of the hard face is observed when combining functional gradients with hierarchical designs, significantly increasing composite properties. Microscopic analysis, digital image correlation, and strain distributions predicted with XFEM are used to discuss the mechanisms responsible for the observed behaviors.
at different scales [18][19][20][21][22] combined with computational models that could reveal the role of different mechanistic ingredients [23][24][25][26] A major challenge in such approaches is the high degrees of interspecies [27][28][29] and intraspecies [30,31] variability as well as location dependency of the measured properties [32][33][34] that make it extremely challenging to perform quantitative and systematic analyses aimed at answering the aforementioned questions and unraveling the mystery of the human bone.
Here, we propose an alternative "voxel-by-voxel" approach, which is analogous to the "atom-by-atom" design paradigm often cited in material science [35][36][37] and aims to unravel the underlying mechanisms of bone through multiscale design of bone-like hard-soft composites. The different "mechanistic ingredients" can then be systematically incorporated in the design of such bone-like composites to study their effects quantitatively and with a minimum degree of variability. Such an approach was impossible to apply until a few years ago because its success is dependent on the availability of multi-material additive manufacturing (AM ¼ 3D printing) techniques that allow for the fabrication of hard-soft composites at multiple length scales. However, recent advances in multi-material AM techniques, including voxel-based polyjet 3D printing, [38][39][40] have enabled a systematic, quantitative, and low-variability study of the roles of different material and mechanistic ingredients in creating the highly efficient properties of hard tissues.

Methods
To separate the effects of different mechanisms from each other, we designed five different experimental groups including 1) a group of specimens with functional gradients (GRAD), 2) a group of specimens with a single-level brick-and-mortar design but without functional gradients (BM-SL), 3) a group of specimens with both a single-level brick-and-mortar design and functional gradients (BM-GRAD-SL), 4) a group of specimens with a two-level hierarchical brick-and-mortar design but without functional gradients (BM-2L), and 5) a group of specimen incorporating both a two-level hierarchical brick-and-mortar design and functional gradients (BM-GRAD-2L) ( Figure 1). The specimens used in all five groups were single-edge notched fracture specimens.
We used a multi-material AM technique working on the basis of jetting multiple streams of UV-curable photopolymers (Objet350 Connex3™, polyjet multi-material 3D printer, Stratasys Ltd., USA) [38] for the fabrication of our specimens. The minimum resolution of the printer was 600 dpi along the x-axis, 300 dpi along the y-axis, and 847 dpi along the z-axis, resulting in a cuboid voxel with the following dimensions: 42 Â 84 Â 30 μm 3 . This cuboid was, therefore, the minimum achievable building block used in our designs. We used commercially available materials VeroCyan (RGD841) and Agilus30 Black (FLX985) for the deposition of the hard and soft phases, Figure 1. A schematic drawing of a single-edge notched fracture specimen a) with the corresponding dimensions used in the current study. We used brick-and-mortar unit cells to create structures with b) one and c) two levels of hierarchy. Binary images were used to 3D print specimens with graded designs, nongraded designs, and a combination of both. d) The specimens with a gradient (GRAD) had 100% hard volume fraction at the regions far from the crack line, which was linearly and symmetrically reduced to 40% in the middle of the specimens along with the initial crack position. e) We used a brick-and-mortar subunit to create our nongraded architected structures (BM-SL). f ) The superposition of the graded and non-graded specimens created a new group (BM-GRAD-SL). g) The second level of hierarchy was also defined by combining several subunit cells (BM-2L). h) The superstition of the GRAD and BM-2L designs resulted in the design of the double-level brick-and-mortar specimens with a gradient (BM-GRAD-2L).  -h, bottom row, and  Table S2, Supporting Information). The target overall hard volume fraction (ρ o ) of the specimens was 71% (Table S2, Supporting Information). The designs were defined through binary images files which acted as the input files for the 3D printer and had values of either 1 or 0 assigned to the hard or soft phases, respectively. To process the designs, we introduced a linear gradient pattern in the specimens with functional gradients and defined it in a gray-scale image. These gray-scale images were later converted into binary images using halftoning algorithms [41] (i.e., GRAD, BM-GRAD-SL, BM-GRAD-2L groups) (Figure 1d,f,h). The specimens without functional gradients (i.e., BM-SL, BM-2L groups) had only brick-and-mortar patterns in their structures, which were directly designed as binary images (Figure 1e,g). The BM-2L specimens had similar brick-and-mortar subunits in their first level of hierarchy similar to those of BM-SL ( Figure 1b). Incorporating those subunits as the building blocks of another brick-andmortar arrangement created a second level of hierarchy for the BM-2L specimens ( Figure 1c). The image analysis steps were performed in Matlab (R2017b, Mathworks, USA).
As there are no standards available to follow for the design of these functionally graded composites, we designed our fracture specimens based on the recommendations provided in the literature for the fracture tests of ductile and brittle composites. [42,43] The geometrical parameters of the specimens are shown in Figure 1. The initial crack was placed in the middle of the specimens, spanned 20% of their widths, and was placed at the weakest point with lowest hard volume fraction. To attach the specimens to the grippers of the mechanical testing machine, the lower and upper sides of the specimens were extended. The extended parts were made from the same hard material as the hard phase of the composites and were printed together with the specimens.
To perform the fracture tests, the specimens were subjected to a quasistatic uniaxial tensile load (stroke rate ¼ 2 mm min À1 ). Unless otherwise stated, a Lloyd machine (LR5K, load cell ¼ 5 kN) was used for the mechanical tests. The specimens were attached to the machine by a custom-made aluminum gripping system and were fixed using aluminum pins. A preload equal to 1 N was applied at the beginning of the tests. The time (T ), force (F), and deformation (u) were recorded by the testing equipment.
The normal stress (σ) was calculated as the applied load (F) divided by the initial effective cross-section area: A 0 ¼ t Â ðW À a 0 Þ, where t is the thickness, W is the width, and a 0 is the length of the initial crack. The strain (ε) was defined as the relative displacement (u) with respect to the initial length of the specimen (L). The stiffness (E) was defined as the stiffest slope of the elastic region of the stress-strain curves. To calculate the stiffness, a moving regression algorithm with a bounding box of 0.2% was used to find the greatest value of the slope. The fracture stress (σ max ) was defined as the maximum stress registered during the tests. The fracture energy (U ) was defined as the energy required to break the specimen and was calculated by numerically integrating the area below the stress-strain curves.
We also measured the full-field strain maps using the digital image correlation (DIC) technique. The DIC tests were conducted on a Zwick Roell machine (Zwick GmbH & Co. KG, load cell ¼ 20 kN). The surface of the specimen was first cleaned and then speckled by black dot patterns. To enhance the contrast with the black speckles, the specimens were illuminated using a profilux LED during the experiments. Two high-resolution digital cameras (4 MP with CMOS chip) recorded the movements of the speckle patterns during our experiments. The full-field strain maps were then determined using the associated commercial software (Vic-3D 8, Correlated Solutions, SC, USA).
We also performed computational simulations of the crack propagation tests using the extended finite element method (XFEM) implemented in a commercially available nonlinear code (Abaqus v.6.14, Dassault Systèmes Simulia, USA). Plane strain quadratic elements (CPE8) were used for the discretization of the specimen geometries defined by the bitmap images, which were also used for the 3D printing of the specimens. Ideally, the geometry of the finite element models has to be discretized up to the size of the voxel size of the 3D printed structures (i.e., 42 μm), as this would allow for matching the manufacturing resolution. However, that would have resulted in a formidably large number of elements for the XFEM simulations. To decrease the number of elements in our computational models, we coarse-grained the bitmap images and assigned an average value of the hard volume fraction to each coarse-grained element. Each element within the XFEM models, therefore, represented a number of neighboring voxels from the original binary image ( Figure 2). The density of each element was, then, calculated using an averaged gray-scale value of all the original voxels it represented. Therefore, the total density of each voxel dictates the final material properties of that element. The new density of each element with the resolution of n bits per edge was calculated as ρ 0 ¼ ρ vox ¼ P bits hard n 2 . This density was used to calculate the individual element's elastic modulus Eðρ 0 Þ, separation strain ϵ sep ðρ 0 Þ for the crack separation analysis. From our experimental data, we obtained polynomial expressions correlating the mechanical properties and the hard volume fraction of the new element (Table S1, Supporting Information). The crack propagation criteria were set based on the maximum strain separation, which means two elements were separated only if their strain values exceeded the maximum strain separation value. These models and their parameters are shown in Figure 2 and Table S1, Supporting Information.
The monolithic experimental data were used to obtain the individual hard and soft material properties (Figure 2d,e), and were also used to obtain the polynomial expressions. We applied displacement control at the top of the model while fully bounding the lower end. The displacement value increased until full separation was achieved. Strain fields and reaction forces were recorded throughout the computational simulations.

Results and Discussion
The stress-strain curves of the specimens in different groups showed a similar trend in which the reaction force (or stress) linearly increased. This trend continued until the ductile fracture phase started (Figure 3a). After reaching their maximum values, all groups exhibited nonlinearly decreasing forces and, thus, softening behaviors. The specimens without gradients (BM-SL and BM-2L) presented the least ductile fracture, as their ultimate strains were the lowest among the different designs (Figure 3a). Contrary to those, the specimens from the GRAD group showed the highest ultimate strain but the lowest maximal stress (Figure 3a,b). Interestingly, the specimens with a combination of brick-and-mortar structures and gradients (i.e., BM-GRAD-SL, BM-GRAD-2L) showed a mixture of the behaviors of both graded and nongraded designs (Figure 3b). The initial crack propagated in a straight line for the specimens with gradients, while crack-bridging was observed for the nongraded specimens (Figure 3c,d,f ). The fracture surface of the GRAD specimens showed a high concentration of microcracks in front of the initial crack ( Figure 3c). As the initial crack propagated in the specimen, these microcracks coalesced and formed curved cracks that deflected to the free surface of the specimens (Figure 3c). However, when the initial crack had almost moved to the middle of the specimen, the majority of the macrocracks were parallel to the fracture surface, exhibiting a flat surface with granular-like profile (Figure 3c).
The fracture surface of the BM-SL specimens showed very different characteristics than those of the GRAD specimens (Figure 3d). Several semielliptical and elliptical macrocracks were observed on the fracture surface of the BM-SL specimens (Figure 3d, SEM subfigures), which could be due to the crack bridging occurred as the specimens ruptured (Figure 3d).
During testing, these specimens failed more swiftly than the others, which correlates with the fact that they were also the most brittle. In addition, the overall surface of these specimens was less flat than those of the GRAD group (Figure 3c,d). Introducing a gradient in the BM-SL specimens resulted in the formation of microcracks that were in parallel with the width of the specimens (Figure 3e). This formation continued as the initial crack propagated in the specimens until it reached the last third of the specimen after which the crack was deflected (Figure 3e). Overall, a considerably more granular profile was observed for the BM-GRAD-SL specimens (Figure 3e), but this profile was not more than the fully gradient designs. These observations suggest that the failure mechanisms of the BM-GRAD-SL group are a combination of both of those mechanisms that were previously discussed for the GRAD and BM-SL specimens. The mortars in the second level of the BM-2L specimens acted as crack barriers preventing further microcrack propagation (Figure 3f, SEM subfigures). Unlike the BM-SL specimens, most of the microcracks that were formed on the fracture surface of the BM-2L specimens were parallel to the crack front (Figure 3f ), and their orientation did not change as the crack propagated. Introducing gradients in these specimens (BM-GRAD-2L), smoothened and diluted the mortar lines observed in the BM-2L specimens (Figure 3g). Therefore, the size of the microcracks in the BM-GRAD-2L specimens was notably smaller than those of the BM-2L group. The microcracks in the BM-GRAD-2L group were mostly distributed in parallel with their widths (Figure 3g). All the failure mechanisms observed in previous designs were also present in the latter group. This synergy of mechanisms may be responsible for the overall higher BM-GRAD-2L performance as compared with those of other design groups.
The mechanical properties of the monolithic and composite materials are compared in Table S3, Supporting Information.  Table S1, Supporting Information. A comparison of the stress-strain curves predicted by our computational models (XFEM) and experimental observations for monolithically d) soft and e) hard materials.
www.advancedsciencenews.com www.aem-journal.com When compared with the structured patterns, GRAD specimens were found to have lower stiffness and maximal stress, but higher fracture strains (Figure 4a-c). However, having an ordered architecture (i.e., brick-and-mortar) in the BM-SL and BM-2L groups could activate other toughening mechanisms in those structures, thereby resulting in higher stiffness and fracture stress values.
Comparing the properties of the groups BM-SL and BM-2L showed that the specimens with a single level of brick-and-mortar hierarchy had slightly higher stiffnesses and fracture stresses than those with two-level hierarchies (Figure 4c). However, the fracture energies of the BM-2L specimens were up to 1.5 and 1.3 times higher than those of the BM-SL and monolithically hard www.advancedsciencenews.com www.aem-journal.com specimens, respectively (Figure 4a,b). This can be due to the activation of the crack-dissipation mechanisms caused by the presence of a thicker mortar line at the second hierarchical level in the BM-2L specimens. In addition, the structured organization of hard/soft materials at the second hierarchical level could be paralleled as an "impurity placement" in the composite, preventing crack propagation. Implementing the gradient in the BM-GRAD-SL and BM-GRAD-2L specimens did not influence the stiffness values and fracture stresses as compared with those without gradients (Figure 4a-c). The fracture energies of the groups with both single-level and two-level hierarchies, however, increased after introducing the gradients (Figure 4a,b). This can be explained by the fact that the brick-and-mortars in the BM-GRAD-SL and BM-GRAD-2L groups did not have purely hard fibers (or soft matrices) in their microstructures. The random distribution of the particles inside the bricks and mortars can create more obstacles along the crack path, which can make it more difficult for the crack to propagate through the specimen. Therefore, more energy was required to break the specimen apart. Interestingly, comparing our results with those of purely hard specimens showed that the fracture energy of the BM-GRAD-2L specimens reached higher values (almost 2 times) than those of the monolithically hard ones (Figure 4a,b).
We also adjusted the mechanical properties of each specimen with respect to their target overall hard volume fraction,ρ o , and www.advancedsciencenews.com www.aem-journal.com the target hard volume fraction in front of the crack tip,ρ c ( Figure S1, Supporting Information). As the difference between the actual and target hard volume fractions were relatively small (i.e., %3.5%, Table S2, Supporting Information), we assumed a linear correlation between the mechanical properties and the hard volume fractions. This adjustment did not change the order of the mechanical properties, and graded structures with two-level hierarchies showed higher fracture energies among all cases ( Figure S1, Supporting Information). We validated the accuracy of our computational models by performing a one-to-one comparison with the average experimental data for every design for different mechanical properties ( Figure S2, Supporting Information). The group rankings between experiments and computational models were consistent for every property regardless of the type of mechanical properties. Our computational models could also predict the values of stiffness and fracture stress with high levels of accuracy. In our computational models, the fracture energies were calculated up to the point that they reached the maximum stress value. Slight deviations from the one-to-one correspondence line ( Figure S2c, Supporting Information) for the fracture energies calculated using our computational models as compared with the experimental values was likely caused by the assumptions we made in our models. Those assumptions involved using 2D-plane strain elements and the coarse-graining of our elements. In addition, we assumed a strain separation value for the propagation of the cracks. Choosing an energy-based criterion for the element separation may, therefore, mitigate some of those deviations.
The strain patterns resulting from the DIC measurements and computational models were very similar (Figure 4d-h). In all groups with and without gradients, a high strain concentration was seen in front of the crack tip (Figure 4d-h). This strain localization in the BM-SL and BM-2L groups formed a butterflyshaped region in the crack tip, characteristic of mode I crack propagation experiments (Figure 4d,f ). The existence of the brick-and-mortar patterns in these specimens did not substantially affect the strain distribution.
The strain distributions for the specimens with a gradient (i.e., the BM-GRAD-SL, BM-GRAD-2L and GRAD groups) showed a more substantial zone of strain localization as well as higher values of the von-Mises strain (Figure 4e,h). The strain distribution for the graded specimens was distributed across the entire crack path. The length of this strain distribution was larger for the graded specimens (i.e., the GRAD group) (Figure 4h). The overall patterns of strain distribution obtained from our computational models were consistent with our experimental observations (Figure 4d-h). These strain distributions confirm that introducing gradients can significantly change the strain distribution in front of the crack tip and directly affect the energy required for crack propagation.
The biomimetic design approaches proposed in this study can have several applications in the field of biomedical engineering such as in the fabrication of functionally graded dental [44,45] and orthopaedic implants [46] as well as engineering complex tissues [47] or interfaces. [48][49][50] The combination of hierarchical designs and functionally graded paradigms can also be used together with optimization approaches to find optimal designs required for better tissue (bone) remodeling. [32,51] We have so far discussed various toughening mechanisms, including granular-like reinforcement (for the graded specimens), multiscale hierarchical arrangements, crack-bridging, and crack barriers. There are also other potential geometry-based toughening mechanisms that were not considered in the current study. For example, introducing certain arrangements of short fibers instead of granular-like reinforcements may add additional toughening mechanisms to the composites. Moreover, different orientations or organizations of bricks-and-mortars unit cells can also be explored in combination with the presented mechanisms. This includes geometries such as chiral, helix, or diamond-like unit cell arrangements. Optimizing the arrangements of such unit cells can result in highly improved properties of advanced architected materials. Finally, in addition to the role of geometry, printing at even finer resolutions (e.g., nanoscale resolutions) using other AM techniques such as direct laser writing (DLW) [52,53] can potentially activate other fracture mechanisms similar to those seen in bone.

Conclusion
In summary, our results suggest that gradients, as a standalone tool, cannot affect the fracture properties of functionally graded hard-soft composites. However, functional gradients with hierarchical designs enable hard-soft brick-and-mortar composites to achieve higher fracture energies as compared with the composites that only benefit from hierarchical designs. These combinations can reach up to 2 times higher fracture energies than the hard phase itself. The results presented here can open up new frontiers in understanding how nature combines mechanistic tools to achieve extremely efficient structural materials and assist in the development of designer materials [54][55][56][57] that have numerous potential applications in high-tech industries.

Supporting Information
Supporting Information is available from the Wiley Online Library or from the author.