Bioinspired rational design of bi-material 3D printed soft-hard interfaces

Durable interfacing of hard and soft materials is a major design challenge caused by the ensuing stress concentrations. In nature, soft-hard interfaces exhibit remarkable mechanical performance, with failures rarely happening at the interface. Here, we mimic the strategies observed in nature to design efficient soft-hard interfaces. We base our geometrical designs on triply periodic minimal surfaces (i.e., Octo, Diamond, and Gyroid), collagen-like triple helices, and randomly distributed particles. A combination of computational simulations and experimental techniques, including uniaxial tensile and quad-lap shear tests, are used to characterize the mechanical performance of the interfaces. Our analyses suggest that smooth interdigitated connections, compliant gradient transitions, and either decreasing or constraining strain concentrations lead to simultaneously strong and tough interfaces. We generate additional interfaces where the abovementioned toughening mechanisms work synergistically to create soft-hard interfaces with strengths approaching the upper achievable limit and enhancing toughness values by 50%, as compared to the control group.

At the same time, we used the simple parametric equations of a helix for the design of the collagen-inspired specimens 2 .We defined these equations in the cartesian system as:  = ;  =  sin( + ) ;  =  cos( + ) … (4)   where  ∈ [0,2) defines the unit-cell discretization,  = 2.032 2 mm is a constant that describes the separation of the helix,  = 0.508mm is the radius of the helix, and  = [0, 1 ]  is the phase angle that defines the rotations of each of the three helices in the design.
After repeating each of these cubic 2.032mm edged unit-cells across the total dimensions of the interface, we used the () function to define the thickness of the helical beams, thereby achieving the proper material discretization.
We processed an additional set of simulations to study the effects of the number of modeled unit cells on the simulation results.In particular, we tried to establish whether the models in which only a single unit cell was included could capture the behavior of the actual experimental tensile tests.To this end, we prepared simulations of 3 different designs.These were the nongraded control (  = 0 mm), the long OC (  = 12 mm), and the long PA (  = 12 mm) (Supplementary Figure 5).We selected these designs because they are quadrant-symmetric, allowing us to simulate only 2 unit cells of the interface (i.e., 1 along their thickness and 2 across their width).These simulations were prepared under the same conditions as described in the main article.However, symmetric boundary conditions were additionally applied on the  = 0 (i.e.,   = 0,   = 0,   = 0 ) and  = 0 planes (i.e.,   = 0,   = 0,   = 0).
The estimations made by the quadrant symmetric (Quad) model of the control group were remarkably similar to those of the single unit cell (UC) models and DIC-measured strain distributions (Supplementary Figure 5a).In both types of simulations, the deformations were highly concentrated at the interface, with the Quad and UC estimations showing   values of 2.28 and 1.7, respectively.Although the absolute strain values differed between both types of models, they occurred at the corners of the specimens regardless of the model type.The corners of the specimens are the locations where high values of shear deformation tended to be present.
These peak strain values also appeared in the DIC measurements, where, as expected, the   value was lower than those predicted computationally.This is due to the averaging effects caused by the lower resolution of the DIC measurements as compared to the computational models and the potential photopolymer mixing prior to curing.Finally, both UC and Quad simulations showed diminished values of the von Mises strains at the center of the cross-sections of the specimens, which is expected because, in this region, only volumetric strains occur.
In the case of the long OC geometries, the strain patterns of the first layer of the Quad and UC simulations were almost identical to those of the DIC measurements (Supplementary Figure 5b).Similarly, the patterns of the  , () distribution within the first layer of the specimens presented the same transitions as observed in the experimental results, indicating that both types of models capture the essential features of the experimental observations.Consistent patterns of strain were also present within the layers with the maximum strain value, regardless of the type of the FEM model.The FEM-predicted  values were 6.36 and 4.53 for the Quad and UC simulations, respectively.Although the magnitude of these strain peaks differed, the shapes of the maximum  , () plots followed the same trends, with the strain peaks occurring in the regions with a hard discontinuous transition from the hard to the soft material.
Similarly, the () functions estimated by both UC and Quad models were similar both in values and trends, indicating that the number of the modeled unit cells does not alter the essence of the predicted mechanical behavior.
Considering the long PA simulations, both Quad and UC models estimated somewhat higher peak strain values within the top layer of the interface than those measured through DIC (Supplementary Figure 5c).As discussed in the main text of the article, the absence of these peaks in the DIC images was likely due to its lower resolution and the potential material mixing between single droplets of the photopolymers.Despite these differences, the  , () values and their localizations were remarkably similar in both the Quad and UC models.Similar to both other designs, the strain values predicted by the Quad models were higher than those resulting from the UC simulations (FEM(max)  = 1.84 and 1.79 for the Quad and UC models, respectively).In both cases, such strains were localized within the middle region of the interfaces and along the edges of the specimens.As described in the main text, the strain peaks observed within the regions made from mostly hard material are of relatively minor importance for the PA designs because any cracks initiated in these regions are likely to be arrested.Moreover, the plots of  , () were remarkably similar in the regions closest to the soft material, indicating that the maximum strains estimated by the UC model are representative of the full specimens at regions where failure is most likely to occur.As in the case of the long OC specimens, little to no differences were present between the trends seen in the plots corresponding to the UC and Quad () models, further corroborating the claim that UC models can successfully capture the main mechanisms behind the observed differential mechanical behaviors of the various considered groups.
Overall, the main differences observed between the UC and Quad models were the FEM(max)  values which were somewhat higher for the quadrant symmetric simulations.
The localization of these peaks was, however, highly consistent between both types of models.
Moreover, the differences in the peak strain values tended to disappear in the regions closest to the soft material, which are the most critical places (because failure due to a sub-optimal interface is likely to occur within those regions).Finally, the similarity between the elastic modulus functions predicted by both types of models allowed us to conclude that UC FEM models are representative of the full-size models, allowing for the application of such computationally efficient simulations for the analysis of the remaining designs.

Supplementary Note 3
An additional set of measurements were performed to compare the three-dimensional behavior of different interface designs under loading and to corroborate the performance of our computational models.For these experiments, we used TESCAN CorTOM micro-CT scanner (TESCAN, Brno, Czech Republic) for micro-CT imaging of a non-graded (Ctrl), a PA (  = 12 mm), and an OC (  = 12 mm) design.To enable in situ assessment of the deformation behavior of the specimens, we scanned the specimens under two loading conditions: without load and under a 6 N load (Supplementary Figure 6a).A load of 6 N was selected because this was the load used for the Quad FEM simulations presented in Supplementary Note 2. The specimens were scanned over a 360˚ rotation with a 3D pixel size of 31 µm and an angular rotation step of 0.25˚.The CT images were acquired at 120 kV and 260 µA with a total duration of ~5 min per imaging cycle.The image analysis of the specimens was performed using the software Fiji (v1.53).Since no internal material pattern was present within the 3D images or any difference was present between the greyscale values of the hard and soft material regions, no digital volume correlation (DVC) analysis could be performed.However, we measured the cross-sectional areas along the length of the specimens and compared with their computational equivalents as per the FEM simulations reported in Supplementary Note 2 (Supplementary Figure 5).
When deformed, the Ctrl design presented a region with sharp deformations at the location of the soft-hard transition (Supplementary Figure 6b).Moreover, its surface area pattern under loading showed an abrupt transition of magnitude at the same region, confirming the presence of high strain concentrations at the edges of the interface.Contrary to these results, the PA and OC designs did not show the aforementioned high deformations at the corners of the interface, and their surface area plots showed a smooth transition at the interface edges.These measured surface area curves were highly similar to those estimated with the FEM simulations with smooth transitions present for the PA and OC designs and the sharp transition characterizing the Ctrl one.Furthermore, when loading the OC design, two gaps were detected in the middle of the interface when loading the specimen (Figures S4d).The locations of these gaps corresponded to the regions where the highest strain concentrations of this study were found, which was at the hard material discontinuity of the long OC interface.Overall, these results confirm that our FEM models and the study of strain concentrations at soft material regions are capable of capturing the internal 3D behavior of the presented 3D soft-hard interface designs.
Moreover, these results call for future studies with in-situ micro-CT scanning under continuous loading, where the 3D-printing equipment is modified to include micro-or nanoparticles that can be visualized with a micro-CT scanner, allowing for the recognition of the internal structure of the specimen and performing DVC.

Supplementary Note 4
We performed additional simulations to determine the failure modes that can occur within various soft-hard interface designs.These simulations were performed for the control (i.e., nongraded), short (i.e.,   = 4 mm) GY, long (i.e.,   = 12 mm) CO, PA, and GP designs with the same meshes as the ones used for the initial hyperelastic simulations.In this case, however, we introduced plasticity and ductile damage with element deletion into the material models used for the soft material, allowing us to study the failure route of each design.These failure considerations were chosen since they yield an efficient and simple method to analyze the failure mode of the interfaces at multiple voxel locations without a need for predefining a single propagating crack or for including cohesive zone elements between every hard and soft voxel of the interface.Due to software (i.e., Abaqus) limitations, such models required changing the elastic material properties of the soft material from hyperelastic to linear elastic (i.e., elastic modulus  soft = 1.3 MPa and  soft = 0.48) and the type of analysis from quasistatic to dynamic explicit.The yield stress point of the soft material was set to 1.2 MPa under a von Mises yield criterion, followed by a linear strain hardening modulus of 0.58 MPa, where element deletion was set for a plastic strain of 5% (Supplementary Figure 10a).These values are based on the existing coarse-graining models for these materials 3 .Regarding the boundary conditions, we applied a displacement to the hard edge of the specimens until mesh separation while applying symmetric boundary conditions to the soft edge of the mesh.
Following the post-processing stage (Supplementary Videos 10 to 14), we found that the stress vs. strain curves of all the designs show similar increases in the stress followed by failure due to element deletion (Supplementary Figure 10b).The long GP design showed the highest values of the strain energy density, while the short GY design underperformed when compared to the other designs (Supplementary Figure 10c).Regarding the plastic energy dissipation prior to failure, the non-graded design showed hardly any capacity for energy dissipation (Supplementary Figure 10d), instead featuring high strain concentrations at the edges of the interface, which caused its sudden failure (Supplementary Figure 10e).Similarly, in the short GY models, the sharp-ended features at the edges of the specimens introduced strain concentrations that caused the separation of the mesh (Supplementary Figure 10e), confirming the design guideline presented in the main text that recommends avoiding such features.It is, however, important to emphasize that the short GY designs showed the lowest strain concentration values among all the TPMS designs in the initial hyperelastic FEM simulations.
Indeed, potential photopolymer blending effects within the 3D-printed specimens may have ameliorated the negative effects of these features, explaining their improved performance in the initial experiments.The long CO and GP designs showed comparatively high strains along their interfaces with regions of relatively low values of strain concentrations that resulted in multiple regions of plastic deformations, confirming how relatively compliant transitions can help in accommodating strains through plastic energy dissipation and can prevent the catastrophic failure of more periodic interfaces.Moreover, both the PA and GP designs showed multiple regions of soft material separation along their interface prior to failure, where cracks were arrested by neighboring hard material voxels.These observations confirm how particlebased designs provide the additional toughening mechanism of crack deflection.However, it is important to mention that for such a mechanism to be applicable, the distribution of the randomly positioned particles must be periodic 4 , unlike in these simulations where particles were present only in a very small cross-section area.All in all, this comparative study allowed us to study how the presented design guidelines that suggest avoiding sharp-ended geometries the zone in which the different photopolymer phases were mixed as well as additional ductile fracture computational simulations that considered photopolymer blending.
The mechanical characterization of the soft-hard mixing zone was performed using atomic force microscopy (AFM, JPK Nanowizard 3, Berlin, Germany) on a non-graded specimen (Supplementary Figure 12a).We utilized a TESP-HAR probe (Bruker, Billerica, Massachusetts, USA) with a length of 125 μm, a width of 40 µm, a thickness of 4 µm, and a nominal spring constant of 42 N/m to measure the elastic modulus change at the transition between the hard and the soft 3D printed material.The calibration of the probe was performed using the thermal tuning method and resulted in a sensitivity of 26.63 nm/V and a stiffness constant of 49.953 N/m.We indented an area of 250 × 100 μm 2 with an indenter separation of 0.39 μm, yielding 250 × 100 indentation force-displacement curves.The indentation curves were acquired for a force setpoint in the range of 1 -4 μN.For each curve, the elastic modulus was calculated using the Hertz-Sneddon model, assuming the AFM tip to be a flat cylinder with a tip radius of 10 nm (the nominal value) with the following relation: , where  is the elastic modulus, is the Poisson's ratio,  is the nominal tip radius, and  is the displacement.After processing the data, a sigmoidal function was used to characterize the mixing zone (Supplementary Figure 12b), similar to what other studies have proposed 5 .This function is defined as: , where () is the elastic modulus along the interface,  ,ℎ and  , are the hard and soft material elastic modulus values, respectively, and  and  are constants. ,ℎ and  , were calculated from the average  value of the farthermost 10 μm from their respective sides of the interface.The values of  and  were obtained by fitting the function using the nonlinear least-squares method.The size of the mixing zone was calculated as the distance between the regions exhibiting elastic moduli corresponding to 0.5% and 99.5% of the elastic modulus of the hard material in the sigmoid design.The obtained size was 151.56 μm, which is in good agreement with the ~150 μm transition reported in a previous study 5 .Moreover, 90% of the material transition (i.e., 5% to 95% hard material) occurs in a span of 84.37 μm, which matches the nominal dimension of each voxel (i.e., 84 μm)..Additional ductile fracture analyses were performed to unravel the potential effects of photopolymer blending on the performance of soft-hard interfaces.Toward this end, we introduced an algorithm that allowed us create an intermediate phase between the soft and hard phases throughout the entire geometry of the specimens in our computational models (Supplementary Figure 12c).For this intermediate phase, we first changed the element-tovoxel ratio from 1 to 27 (3×3×3).To consider the potential blending at the edges of each voxel 5 .We also updated the value of the ratio of the hard material '  of each element across the ijk directions (i.e., initially  ′ = 0 for pure soft and ' = 1 for pure hard material) by averaging it with its adjacent elements.This averaging was performed using the following equation: The resulting ' values of each element, which were in accordance with the mixing zone fit presented above (Supplementary Figure 12d), were then used as input for obtaining the mechanical properties of each element.For this purpose, coarse-graining models for large deformations similar to those existing in the literature were used as the constitutive models 3 .
We applied this process for the non-graded control design and the long PA ones and compared their resulting mechanical response with those of their non-blended variants.Since increasing the element to voxel representation by 27 would make the computational models infeasible to run, we restricted the meshes to interfaces of 12×12 voxels in the surface area and 96 voxels in length, yielding 373248 C3D8 elements per mesh.The resulting four simulations were processed under the same conditions as described in the Supplementary Note 4 of this document.
After post-processing the results (Supplementary Videos 15-18), the blended version of both the control and the PA gradients showed improvements in their ultimate stress and strain energy density (Supplementary Figure 12 e-g).For the non-graded design, the strength increased from 0.382 MPa to 0.399 MPa, while the strain energy density increased from 20.04 kJ/m 3 to 22.37 kJ/m 3 .As for the PA, the strength increased from 0.41 to 0.45 MPa, whereas the strain energy density improved from 35.96 to 44.6 kJ/m 3 .The relatively larger improvements observed in the particle-based design can be partially attributed to the increased contact surface area between the polymer phases.It is, however, important to acknowledge that the separation of the blend phases highly depends on the selected coarse-graining model.In this case, we assumed that the ultimate strain before the separation of the intermediate phase is given by the rule of mixtures between the hard and soft material phases.This assumption was necessary because no experimental data exist regarding the ultimate strain of the sub-voxel blended phase.That said, the overall strain distribution of the blended simulations and the failure modes of both designs were remarkably similar regardless of whether blending was implemented in the models (Figures S13h).In the case of the PA designs, multiple non-critical cracks were present prior to failure for both simulations, indicating that the crack deflection capacity of this design remains present regardless of whether blending is implemented in the models.
In short, these analyses showed that while photopolymer mixing plays a potential toughening role by reducing interfacial strain concentrations, the overall behavior and the toughening mechanisms achieved by varying the geometry of the interface remain consistent between the models incorporating such a blending effect and those neglecting it.To properly account for the additional toughening mechanisms resulting from blending, one would need to perform additional experimental characterizations of the effects of photopolymer blending on the complex particle-based composites, building up on this and other existing analyses of soft-hard slab connections 5,6 .
the length of the gradient as well as permuting the interface geometry present at the end of the crack tip.
We defined a symmetrical hard-soft-hard interface geometry for the fracture toughness test specimens, with nominal dimensions of 24.4×24.4mm 2 in length and width (L×W), as well as an out-of-plane thickness of 24.4 mm (B) to promote critical crack propagation (Supplementary Figure 13a).The crack, with an initial length of  0 = 6.1 mm, was prescribed at the soft end of the interface, which was the location where cracks propagated for most of the inefficient interfaces under uniaxial tensile deformations.The selected FG geometries were GY, PA, GP (  = 8.13 mm), and a non-graded control (Ctrl).After 3D-printing three specimens for each design, we tested them for mode I deformations using the same equipment and setup as the tensile tests described in Section 2.3 of the main manuscript.For post-processing, we calculated the fracture toughness of each test in terms of the J-integral,  0 , utilizing the standard for ductile homogeneous materials ISO 12135.This standard defines  0 with the following expression: , where   is the maximum recorded force of the test,  is the elastic modulus (i.e., the slope of the stress-strain curve of the crosshead measured between 0 to 50% maximum stress), and  is the Poisson's ratio of the specimen which is assumed to be 0.45, for simplicity. 2 and   are geometrical factors defined by the following expressions: The remaining term,   , represents the plastic strain energy of the system and is defined in terms of the total strain energy,  , and the elastic strain energy,   , with the following expressions: , where CMOD is the crack mouth opening displacement of the specimen, measured with digital extensometers within the DIC software,  is the slope of the F-CMOD curve measured between 0% and 50% of   , and U was measured as the area under the F-CMOD curve.
After post-processing the data, the results showed a distinct behavior between the non-graded specimens and the FG designs (Supplementary Figure 13b).The non-graded design showed the lowest mean  0 values of the study (i.e., mean = 64.07,SD = 7.6 KJ/m 2 ).These designs were followed by the GP (i.e., mean = 70.22,SD = 5.9 KJ/m 2 ), the GY (i.e., mean = 73.79,SD 3.9 KJ/m 2 ), and the PA graded design (i.e., mean = 78.91,SD 7.5 KJ/m 2 ), which showed the highest performance.These indicate that introducing a graded interface improved the performance of the specimens by 9.6-23.1%.However, it is essential to address that, when inspecting the equivalent von Mises strain fields of the interfaces (Supplementary Figure 13c), only the control group specimens presented a secondary crack that initiated at the interface edge that was not included in the original design.This observation corroborates the observation that the unwanted strain concentrations at the edges of the soft-hard connections play a critical role in the performance of soft-hard connections.Overall, this analysis contributes to the definition a suitable geometry for studying the fracture behavior of soft-hard interfaces and confirms the improved behavior of our transition geometries as compared to their gradient-less equivalents.However, future studies on the configuration of the fracture test specimen require testing configurations that guarantee no additional cracks are propagated during the analysis.
Since it is important to prescribe the crack at a position within the interface where the outcome performance is expected to be minimal, performing a permutation analysis on the crack position through computational techniques may be necessary.We performed a pilot analysis with such a method using linear elastic XFEM simulations.These XFEM analyses comprised two different sets of tests.For the first one, we varied the position of the prescribed crack along the interface and estimated the corresponding stress intensity factor,   (Supplementary Figure 13d).For this purpose, the crack position was permuted between 5 voxels behind the interface (i.e., the soft material side) to 6 voxels within the interface geometry (i.e., toward the hard material side).For the second set of tests, the crack position was fixed at the edge of the interface geometry while the geometry of the FG was shifted between 0 and 12 voxels along the direction of the crack (Supplementary Figure 13e).Such a shift in the interface geometry allows us to evaluate how much the local material configuration ahead of the crack tip affects the   value.The selected interface designs for testing included OC, DI, GY, CO, PA, and GP, with the addition of a non-graded specimen as a control group.This yielded 84 XFEM simulations for the first set of tests and 91 simulations for the second.To make the computational models feasible, we discretized only one interface quadrant (i.e., 144×144 voxels) and only one unit cell in thickness (i.e., 24 voxels) of the fracture toughness specimen design, leading to 497664 hexahedral elements with full integration per simulation.Symmetric boundary conditions were applied at the soft material end (i.e.,   =   =   = 0), at the edge opposite to the interface (i.e.,   =   =   = 0), and at the back of the specimen (i.e.,   =   =   = 0), while a linear stress of 0.186 MPa was prescribed at the hard material end.
After post-processing the simulations, high variations of   were present when changing the crack position along the interface length (Supplementary Figure 13f).The maximum   value of the non-graded design was 1086 MPa.mm 1/2 , estimated at one voxel within the hard material region, which was similar to the value at the edge of the interface, which was 924 MPa.mm 1/2 .
Contrary to this, the maximum estimated values of the FG designs varied between 472 and 646 MPa.mm 1/2 , with most of these estimations being found for cracks located well within the interface geometry.Moreover, when positioned within the interface region, some of the   values of the FGs were as low as 3.76 MPa.mm 1/2 .These differences of up to two orders of magnitude indicate a high dependence of fracture toughness value on the positioning of the crack.Similarly, shifting the FG geometry while leaving the crack stationary showed substantial variations of   (Supplementary Figure 13g).For this case, the behavior of all the geometries was within the range of 3 to 6 MPa.mm 1/2 .However, one FG design can have a lower or higher   than the other designs depending on what part of its geometry is directly ahead of the crack tip.For example, a lattice shift of half a unit cell (i.e., 12 voxels) of a CO design can make its   increase from 3.8 to 5.27 MPa.mm 1/2 and change its rank from the most efficient geometry to the fourth most efficient.Therefore, these results corroborate our hypothesis that measuring the most critical   of any interface design would require extensive permutation analyses to test all the possible configurations of crack positioning within a geometry, making this a highly laborious evaluation method.Consequently and as suggested in the main text, analyzing the performance of FG designs in terms of the strain concentration factor is a far more efficient methodology than the one presented in this section.

Table 1 .
The quasi-static tensile test results corresponding to the various types of functional gradient (mean ± standard deviation).

Table 2 .
The quasi-static shear test results corresponding to the various types of functional gradients (mean ± standard deviation).