Interlaced Topologically Interlocking Lattice for continuous dual-material extrusion

We propose a simple and effective 3D lattice consisting of interlaced horizontal beams in vertically alternating directions which interlock topologically: the interlaced topologically interlocking lattice (ITIL). It


Introduction
Multi-material extrusion 3D printers unlock a plethora of applications through combining the unique material properties of various materials.However, depending on the combination of materials the adhesion between the materials can be excessively weak.While only a small number of all possible 3D printing material combinations may exhibit any incompatibility issues, it is often precisely those incompatible combinations where the different chemical make-up produces interesting applications.For example, polypropylene (PP) is semi-flexible and fatigue-resistant, but has a very weak chemical bond to typical rigid filaments such as polylactic acid (PLA).PP is often used for living hinges, which consist of a single part rather than moving components.Combining these materials unlocks compliant mechanism applications such as visualized in Fig. 1.In such cases it is necessary to rely on mechanical interlocking to prevent the materials from breaking apart from each other.Dovetail interlocking is a common strategy to affix two bodies together; one example can be found in jigsaw puzzles.The pieces interlock and stay connected because of the material stiffness.However, the interlock can be broken in plane by deformation of the pieces, and if the dovetail is not accompanied by horizontal features above an below the pieces could be disassembled by translation alone.See Fig. 2(a).Moreover, because the dovetails widen towards their tip, they cannot easily be printed with continuous extrusion toolpaths.
Therefore, we propose topological interlocking , which is a type of interlocking where the interlock is preserved under continuous deformations such as stretching, twisting or bending of any magnitudethat is: the bodies can only be unlocked by discontinuous deformations such as fracture.In order to achieve topological interlocking both materials have a high genus topology: the holes or tunnels in one material are filled with the other material and vice versa, similar to how the rings of chain mail are linked together.A topological interlock remains effective under any deformation of the base materials and can only be broken by failure in either material.This interlocking principle is robust especially when flexible and deformable materials are concerned.
However, most topologically interlocking geometries would introduce discontinuities in the extrusion process when sliced into layers https://doi.org/10.1016/j.addma.2021.102495Received 7 September 2021; Received in revised form 29 October 2021; Accepted 14 November 2021  The dovetail can be disassembled by translation; the cubic lattice causes highly discontinuous extrusion on some layers; ITIL solves these problems.
for 3D printing, as each slice would contain disconnected islands.See Fig. 2(b).Such discontinuities can easily lead to defects, which influence the dimensional accuracy and the mechanical properties of the resulting part.We therefore have to generate interlocking geometry for which the layers consist of long continuously connected areas for both materials.
It seems impossible to generate topologically interlocking geometry with holes or tunnels while enforcing continuous extrusion for both materials; if the one material leaves a hole in a layer then filling that hole with the other material will cause it to be disconnected from the other regions of the second material.The ITIL lattice consists of long horizontal beams which ensure continuous extrusion, which alternate in orientation along the Z axis in order to connect all beams together, thereby creating linked hoops which constitute the topological interlock.Where the beams of the two orientations meet they form long vertical pillars.However, the small island cross-sections of these pillars are all located in between consecutive layers, so that they never result in disconnected islands.See Fig. 2(c).

Related work
Multi-material additive manufacturing.Multi-material additive manufacturing has unlocked a plethora of applications, by making use of functionally graded material properties, tailored composite materials or multi-material designs.Several review papers cover a wide range of techniques on these topics [1,2].Using multiple materials to create colored surface imagery is commonly performed using MultiJet technologies, but can alternatively be performed using inkjet techniques [3] or even DLP resin printers [4].Such techniques can even be extended to deal with translucency [5] and gloss [6].Several color printing techniques have also been proposed for MEX [7][8][9].
Besides visual attributes, multi-material systems can be used to create parts with graded material properties, by generating composite structures with varying densities of soft and hard materials [10].Finetuning the small-scale geometry in which such materials are deposited can give a more sophisticated control over the induced material properties, [11,12] and adjusting the small-scale geometry throughout the product on a meso scale can increase the performance of the product even more [13].
Some MEX systems for multi-material extrusion have been proposed which operate by extruding multiple materials out of a single nozzle, e.g. by routing multiple filaments into a single mixing nozzle [14] or by creating a single strand of multi-material filament [15,16], but such systems can exhibit hardware issues when the different materials require vastly different processing parameters.Therefore the more upscale multi-material extrusion systems use a separate nozzle for each material [17].
Lattices.Lattices such as beam lattices, triply periodic surfaces and ordered dithering have widely been studied for their mechanical properties.Several review papers provide a comprehensive overview of lattices, which are also known as 'microstructures', 'meso-scale structures', 'cellular materials', etc. [18][19][20] Single material lattices can constitute auxetic materials, or light weight structures with tailored material properties, through fine adjustments of the lattice geometry.Multi-material lattices inherently exhibit topological interlocking, which makes them good candidates for interlocking [21].However, optimizing lattices for adhesion between incompatible materials while adhering to MEX manufacturing constraints has been studied scarcely.
Adhesion.Important factors for adhesion between polymers are entanglement and dissipation [22].The adhesion between layers produced within a body of a single material produced by MEX can be influenced by various process parameters [23,24], as well as the geometry of the toolpaths [25][26][27].The adhesion by which two bodies of different MEX printed materials stick together can be influences by a wide spectrum of pre-treatment methods, process parameters and material properties [21,28].Increasing surface roughness might improve the adhesion between materials [29,30], but this supposed benefit is contested [22].For MEX one could try mixing the materials by overlapping their toolpaths to increase adhesion or create simple straight protrusions in order to increase the friction between the two materials [31].However, if such protrusions bulge outward the adhesion does not merely increase because of the increased friction, but also because it would constitute a dovetail type of interlocking structure.
Interlocking.An interesting type of interlocking can be found in interlocking assemblies where ''locking of an individual element within the assemblage is furnished by the kinematic constraints provided by its neighbours by virtue of the element geometry and the mutual arrangements of the elements within the structure'' [32].Although this concept has also been referred to as 'topological interlocking' in literature, it is rather different from the type of topological interlocking proposed in this paper.This paper pertains to bodies of different materials, and the interlocking between these materials is preserved under continuous deformations, i.e., topologically.While interlocking assemblies lock orthogonal movement when the elements are globally constrained transversely, the interlock is nullified when no such global constraint is present.This limits the application area of interlocking assemblies.
Textiles exhibit a different type of interlocking, where individual strands are woven together into what can be viewed as a lattice of knots.Knot theory is a vast area of research within the mathematical discipline of topology.The mechanical properties of textiles are also thoroughly studied [33].The concept of weaving can be exploited in MEX: weaving extruded strands together across several layers can improve layer adhesion within a part of a single material [34].However, because the strands are extruded from the top, no topological knots can be achieved in MEX.It is not possible for MEX to weave strands into an interlocking textile lattice.
Literature on interlocking patterns for adhesion between incompatible materials in MEX 3D printing has focused on extended 2D interlocking dovetail type of designs, such as jigsaw shapes [35], trapezoidal sutures [36], T-shapes [37,38] and star shapes [39] in the horizontal direction as well as in the vertical direction across several layers [27].Topology optimization can generate complex 2D dovetail interlocking shapes, which fit to the specifics of the design locally [40].However, the dovetail shapes are often relatively large, which limits their applicability and if they are shrunk the widening of the dovetails would introduce discontinuous toolpaths.If such interlocking designs are considered on their own without horizontal features above and below they would allow for disassembly in the vertical direction.By considering aspects of the interlock along the Z axis, the dovetail interlocking concept can be expanded into 3D interlocking structures.
The jigsaw idea can be expanded into a 3D interlocking structure by protruding not only sideways, but also vertically, resulting in a shape resembling a tree [41].Similarly the T-shaped interlocking design with horizontal bars could be expanded with bars in the vertical directions.However, such designs violate the semi-continuous extrusion requirement because the layers above and below the base of the  would contain separated islands of one material, which are difficult to print accurately when the two materials do not adhere to each other.
These issues can be addressed by generating a repeating lattice where the interlocks are connected together; simple straight I shaped extrusions can be linked together by cross beams in order to form a topologically interlocking structure resembling a ladder [42].Such a topologically interlocking structure satisfies the continuity constraint only for a single material and is therefore interesting for applications where the one material is 3D printed, while the other is overmolded silicone.However, a different lattice is required if the extrusion continuity constraint is to be met by both materials.

Contributions
We propose an interlocking lattice consisting of interlaced horizontal beams which satisfy the extrusion continuity constraint.The lattice interlocks all degrees of freedom in 3D space and as such lends itself to interfaces between two bodies of different material of arbitrary orientation.We consider the situation of a horizontally applied tensile force to bodies with a vertical interface, and optimize the lattice in two different orientations w.r.t. the interface using a simple analytical model.The optimized structures are validated using numerical simulation as well as physical experiments.Results show that the ITIL lattice on a vertical interface performs comparably to dovetail type interlocking geometries.Optimizing the lattice for interfaces of arbitrary orientation is left as future work.

Method
This paper considers topologically interlocking structures satisfying extrusion continuity constraints, which will be optimized for tensile strength; the geometry of the designs are optimized to yield or break at the highest horizontally applied tensile stress.The highest administrable force applied to a unit cell of any interlocking structure will be carried by both materials  and .In the plane where tensile failure of material  occurs the total force  will be distributed over the cross-sectional area   ; likewise for material .The theoretically optimal interlocking lattice would be such that the force  will be homogeneously distributed over   and   in such manner that both materials will fail at the same ultimate force:  =  y,   =  y,   .The ultimate force before material  breaks can be improved by increasing the area   , but this is at the cost of material  (and vice versa), since their combined area is limited to the total cross sectional area  total of the lattice.The ultimate strength of any multi-material lattice is therefore limited to: This formula for the ultimate tensile strength of interlocking lattices is a theoretical upper bound.In practice the worst cross-sectional area   is in a different plane than the worst cross-sectional area   , so they do not add up to  total .For example, in the dovetail geometry (Fig. 2(a)) the highest stress of a dovetail of material  is at its base, while the highest stress in  will be at the end of the dovetail of material .The geometry by which the interlock is secured reduces the strength compared to the theoretical upper bound of the ultimate tensile strength.The upper bound can therefore never be reached.

Interlaced topologically interlocking lattice
The topologically interlocking structure we propose consists of horizontal beams alternating in material.On top of that we print another set of beams rotated about Z and on top beams of the first rotation again, etc.Long horizontal beams assure continuous extrusion, while the alternating direction of the beams assures the interlock.See Fig. 2(c).
When joining bodies of different materials, some region around the interface between the two bodies should be replaced by a number of unit cells from such a lattice structure.The length of the transition region  is limited by the design of the product within which the lattice is to be employed.We consider a design constraint of  max = 12 min, = 3.6 mm; this is enough to have some design freedom for the lattice, while limiting the impact on the rest of the product.
We consider a single layer of ITIL cells along the interface.Adding more cells in the direction orthogonal to the interface would make them protrude farther into the two bodies, which limits design applicability.If the two layers of cells would have the same geometry then we would not expect any gain in the ultimate strength of the interlock.Optimizing the geometry of the two layers of cells separately could lead to improved ultimate strength, but falls beyond the scope of this paper.
According to the reasoning above we should fill the interface with a single layer of cells, but how should those cells be oriented w.r.t. the interface?While rotating the lattice about X or Y is impossible due to the continuous extrusion constraint, we are free to rotate about Z.However, given that the tensile stress is applied in a single direction it only makes sense to consider two orientations: straight and diagonal.
For the materials chosen by this study, Ultimaker Tough Polylactic Acid (PLA) and Ultimaker PolyPropylene (PP), the theoretical upper bound comes out to be 8.6 MPa.The structure may not only be subject to tensile failure, but also to shear failure modes.Because of the layer-wise build-up employed by MEX, the tensile properties in the Z direction are different from those in the horizontal plane.See Table 1.The structure is furthermore subject to manufacturing constraints determined by the nozzle sizes and the layer thickness.The standard layer thickness ℎ min is 0.1 mm and the minimal line width  min for a 0.4 mm nozzle is 0.3 mm.
If all dimensions of the beams are set to minimal and the angle between the beams is set to 90 • , the finger of the straight ITIL variant takes up a quarter of the area at the base of the cell, so we can expect the strength to be close to 1 ∕4 y, ≈ 2.6 MPa.For the diagonal ITIL variant the beams take up half of the area at the base of the cell, which would be close to 1 ∕2 y, ≈ 5.3 MPa.However, these ultimate strengths can be improved upon considerably by optimizing the geometry.This section considers these two ITIL variants and analytical models to optimize them.

Straight ITIL variant
A unit cell of the straight ITIL variant is visualized in Fig. 3(a).A cell consists of a single finger of height ℎ f protruding from the body of material  outward and part of a cross beam of height ℎ c angled at 90 • .In order to print the relatively short fingers using continuous extrusion, we employ the constraints that   ≥ 2 min, (where  is either material  or ), so that the toolpaths for the outline of each layer can go back and forth along the finger without interruption.The cross beams are long and continuous enough, so they could be printed using a single extrusion path:   ≥  min, .We set the minimum height of the beams to twice the layer height, so as to be able to recover from manufacturing inaccuracies 1 : ℎ f ≥ 2ℎ min and ℎ c ≥ 2ℎ min .

Straight ITIL variant (whole)
In order to optimize the straight ITIL variant for a maximal tensile strength we consider three types of stress, related to three types of failure mode for either material : tensile stress  XX, for   XX , cross beam shear stress  XZ, for    XZ and Z shear stress  XY, for   XY .See Fig. 3(b).These three types of failure mode for either material are modeled using classical beam theory as having a homogeneous stress distribution.Because the total force  is modeled as being homogeneously distributed over the whole cross beam, the two shear stresses obtain only a portion of the total force.Furthermore, they are divided by two because the beams and pillars are fixed on both sides.See Eqs. ( 2) to (4). 1 If the structure would consist of alternating geometry each layer, then the inaccuracy of the one layer can cause over-extrusion in the next layer, which snowballs the problem upward during printing.

Stresses in the straight ITIL variant
where  ∈ {, } and ¬ =  and ¬ =  Model of the straight ITIL variant subject to where Eqs. ( 10) and ( 12) are duplicated for both materials  ∈ {, } The tensile stress of the whole cell is given by  ∕(  +  )(ℎ f +ℎ c ). Combining the above we obtain the constrained optimization problem given by Eqs. ( 5) and (12).The √ 3 in Eqs. ( 11) and ( 12) comes from the von Mises yield criterion.Although one might consider combining the stresses together using the same criterion, this does not increase the accuracy of the model, because the failure mode planes do not overlap.
Because the objective and all stresses are invariant under various scaling operations, we can choose the value of a subset of the design variables.Because of invariance to uniform scaling we employ the design constraint; scaling up the width and force only increases the cross beam shear values  XZ, , so we employ the minimum width constraint; scaling up the height and force only increases the Z shear values  XY, , so we employ the minimum height constraint.This way we limit the design space from six to three dimensions, i.e.   , ℎ f ,   : Using the formulae of the constrained optimization problem one can find the maximum force and thus the maximum stress of any design; we can rewrite each of the mechanical constraints Eqs.(10) to (12) to give a formula for  and the lowest value of those formulae will give the active failure mode for a given design.The resulting response surfaces are shown in Fig. 5(a).The optimum is at the intersection of four constraint surfaces, which is remarkable for a 3D space.

Broken cross beams model
A careful analysis of the geometry will show that if shear failure   XZ has occurred, there is still interlocking between the two materials.If a part of the cross beams of  has sheared off, still the pillar of material  remains, which is surrounded by material .See Fig. 3(c).Once failure mode   XZ has occurred, still any other failure mode has to occur for the interlock to fail completely.Since the failure mode can only happen by part of the PP cross beam pushing against the part of the PLA cross beam which is being sheared off, both cross beams shear together.Because PLA breaks at a lower strain than PP (see Table 1), we know that for these materials   XZ never occurs unless   XZ has occurred.We will construct a separate model for analyzing the case where the PLA cross beam is broken into segments.
Because part of the cross beam is missing the stress on the remaining part increases.Moreover, the cross beam shear constraint for PP ( cb ) comes into play and the Z shear constraint  zb and the cross beam shear constraint  ca are dropped because of the change in the cross beam of .The shear constraints Eqs.(11) to (12) are replaced by Eqs. ( 14) and (15).

Straight ITIL variant -shear constraints for broken case
The response surface of this model is shown in Fig. 5(b).In order to estimate the ultimate force of a given design one must consider both of these models: the one for the whole and the one for the broken situation.When considering designs where the cross beam constraint  ca is active in the whole model, the maximal force applicable is the highest of the two maximal forces according to the two models.

Diagonal ITIL variant
Besides the straight variant of the ITIL lattice, we also consider the variant where the beams are oriented diagonally to the interface surface.Fig. 4(a) shows a simple cell of the ITIL lattice in the diagonal orientation.Because the stress applied is homogeneous and precisely normal to the interface the optimal structure should be symmetric.Due to symmetry there is no distinction between fingers and cross beams to be made in this model.Both beams should have the same height 2ℎ min Fig. 5. Maximum strength according to analytical models for the straight ITIL variant along three 2D slices of the 3D design space.Four constraint surfaces meet at the optimum.By taking the maximum values of these two models the maximum stress before separation can be calculated.and the angle  between the beams and the interface of both beams is the same.The remaining design variables are:   ,   and .
The simple cell consists of protruding fingers and triangular dents (see red triangles in Fig. 4(a)); however, the advantage these dents give to the strength of the pattern is vastly outweighed by their effect on the total length .We therefore trim the triangular ends on the pattern to arrive at the model shown in Fig. 4(b).However, doing this introduces some sharp edges, with a diameter below the minimum feature size  min .We therefore round the vertical edges using a radius  = 0.15 mm and define the dimensions of the model such that the rounded edges at the ends of both fingers are aligned.See Fig. 6(a).
Because of the diagonal geometry we expect that the stress distribution throughout the structure will be quite complex, and the failure will depend on the deformation of both materials during stretching.Nevertheless, we provide a simplified analytical model assuming the stress is homogeneously distributed and disregarding the influence of deformation on stress distribution.We decompose the force  into one component parallel and another orthogonal to the direction of the beams.We then use the parallel component to determine the tensile stress and the orthogonal component to determine the shear stress in the beam, which combine into a single constraint using the von Mises yield criterion in Eq. (20).The diagonal ITIL model is then defined by the constrained optimization problem given by Eqs. ( 16) to (21).

Model for the diagonal ITIL variant
subject to wb ∶  ≥ 2 min, ( 18) where where Eqs. ( 20) and ( 21) are duplicated for both materials  ∈ {, } Again the stresses are all scale-invariant, so we set   = 2 min, .However, because changes in   given the same  will change the angle  of the beams, we cannot assume the design constraint  d is active.
The resulting response surface can be viewed in Fig. 6(b).Because the height of the fingers is minimal the Z shear constraints are both dominated.

Results
In order to validate our analytical models we compare its predictions against both simulation results and physical tensile tests.While the physical tests constitute the final arbiter on the matter, our resources for performing physical tests are limited.Simulations can easily be performed by running a script for multiple days without user interaction.Given that the simulations make use of the same material properties which were acquired from tensile tests performed by Ultimaker, the simulations can teach us about the validity of the homogeneity assumptions in the analytical models.Moreover, the physical test results are afflicted with a spread in manufacturing inaccuracies.The simulations can therefore enrich the understanding we gain from physical experiments.

Numerical simulation
In order to simulate interlocking structures with a range of design parameters we automatically generate INP files in Abaqus CAE (2020) using a script.Solving an INP file gives us the force-displacement graph, from which we can determine the ultimate tensile strength of a particular design.In order to simulate accurately we used the stressstrain curves from tensile tests on the base materials printed flat on the build platform as the plasticity in tabular form.The simulations were performed in the Abaqus/Explicit solver where the Dynamic, Explicit procedure step was used with a mass scaling factor of 10 7 , using geometric non-linearity and general contact (explicit) to disregard friction for simplicity.
The repeating nature of the interlocking patterns was captured by modeling half of the unit cell and applying symmetry constraints to the  A grid search was used to measure the influence on the ultimate strength along each of the design variables   ,   , ℎ c , along with the total length .The search space was therefore 4D and 2D for the straight and diagonal ITIL model respectively.In order to estimate the optimum we fit a smooth response surface to these data points using a radial basis function (RBF) network [43], with a smoothness of  = 1.

Straight ITIL variant
In order to prevent stress concentrations and adhere to manufacturing accuracy, the vertical edges of the straight ITIL variant were rounded with  = 0.15 mm; see Fig. 7.
Newton's method was used to determine the optimum, starting from the best sampled point.This step only considered the dimensions   and   , because  max is given and ℎ f has to be an integer multiple of ℎ min .The resulting hypersurfaces are visualized in Fig. 8.The obtained optima are shown in Table 2.
We compare these results to our analytical model by adjusting the analytical model to capture the inaccurate Z strength used in the simulations:  yZ, ∶=  y, .See Fig. 9.We then observe that our analytical model on average predicts only 7.8 % higher ultimate strength values than then the FEM simulations, with a standard deviation of 16.2 %.

Diagonal ITIL variant
Modeling the diagonal ITIL variant in Abaqus can be quite cumbersome, since it does not natively support periodic boundary constraints.Whereas this problem can be overcame in the straight ITIL variant because it is symmetric, the diagonal variant is only rotationally symmetric.While a symmetry constraint can be used on the top and bottom, the two sides of the design are mirror images of each other, but also flipped vertically.
However, since the height of the beams is relatively low compared to their width we have observed that the stresses and strains are quite similar in the top and bottom.If we model half of the diagonal ITIL cell by cutting it vertically and apply symmetry constraints to the sides, the induced error is only ±10% compared to simulating an interface consisting of two whole cells.
The results of these simulations are shows in Fig. 10(a).We compare the simulations against our analytical model in Fig. 10(b).The analytical model predicts 0.4% lower ultimate strength values on average with a standard deviation of 10%.

Physical experimental tests
Tensile tests were performed on an Instron 3366 Universal Testing machine at 5 mm∕min.Prints were manufactured on a Ultimaker S5 systems in 5-fold with Ultimaker Green Tough PLA and Ultimaker PP using the default 0.1 mm layer thickness profile, with 100 % infill and a custom brim to make sure both materials stick to the build platform.For PP we print the outer before the inner walls so as to improve the dimensional accuracy.In order to deal with the various widths of the beams we generate toolpaths from STL 3D models using the Cura Arachne Engine beta release [44], which implements a framework for generating variable line width toolpaths to fill small geometry of arbitrary dimensions [45].The Inward Distributed and the Distributed strategy were used on PLA and PP respectively.See Fig. 11.

Model parameters
Straight ITIL.The straight ITIL variant suffers from the curse of dimensionality; even when setting   = 0.6 mm, ℎ c = 0.1 mm and  = 3.6 mm, there are still the three free design variables   ,   and ℎ f to determine.With 5 specimens per sample point and limited resources, the total number of data points we are able to test is limited.We therefore chose to sample close to the two optima of the analytical models: whole and broken, as well as deviations from those optima in both directions along the axis of each design variable of 0.3 mm in   and   and 0.2 mm in ℎ f .See Fig. 12(c).
Each sample of the straight ITIL variant has 5 × 5 cells.Because the repetition of cells is broken at the sides of the specimen, the boundary cells are adjusted for manufacturability and stability.The specimens end with a PLA finger on both sides and in cross beams on both top and bottom.See Fig. 12(a).Dovetail interlocking.We compared our interlocking structures against two interlocking designs: trapezoidal sutures and jigsaw interlocking.See Fig. 13.We used   = 2 min, ,  = 0.3 mm and  = 2.4 for the trapezoidal suture, and   = 2 min, and  = 35 • for the jigsaw interlocking design.We printed samples with both   = 3  and   =  y,∕ y,   = 4.48  .We used 6 and 4 repetition respectively and a height of 5 mm.
Note that the jigsaw interlocking structure is quite similar to the trapezoidal suture, with the addition of semicircles to the ends of the trapezoids.The total length  of the jigsaw structure is 2.96 mm and 4.05 mm for the two   values respectively, so the  max constraint is violated by that structure.
The boundaries of these two structures end in half a PLA lobe, because that is the stiffer material.In order to meet the minimum width constraint there, the sides of the specimen are extruded by  min, .However, because the dovetails wider towards the tip, they require an extra toolpath to be filled densely, which is disconnected from the other toolpaths, thereby violating extrusion continuity constraints.See Figs.11(c) and 11(d).

Physical test results
After tensile testing we can observe various failure modes, such as in Fig. 16.The tensile tests performed result in force-displacement graphs, from which the ultimate tensile strength values are derived.The slope of these graphs after the optimum has been reached furthermore tells us something about the failure mode by which the sample has failed.See Fig. 15.
Because the boundary cells deviate from the regular pattern, computing the ultimate tensile strength can be done in two ways, giving rise to two statistics.We compute the cell stress by dividing the maximum force of the force-displacement graphs by the number of cells and then divide it by the cross sectional area of the cell.We compute the total stress by dividing the force by the total cross-sectional area of the sample, including the extra geometry at the boundaries of the sample.Ideally we would compensate for manufacturing inaccuracies by using measured dimensions of the specimens, but measuring the internal geometry of the interlocking structure is practically infeasible, so we use the dimensions of the 3D mesh instead.The results from the physical experiments, along with the analytical and simulated predictions are gathered in Fig. 14.

Validity of analytical models
When comparing the tensile test results to the analytical models and the simulation results (Fig. 14), we observe that the whole straight ITIL variant is predicted rather well by the models.The analytical model tends to overestimate the strength because of homogeneous stress assumptions, while the simulations tend to underestimate the strength.Although the standard deviation is too high to capture the differences in response with significance, the analytical model follows the trends of the FEM simulations rather well and the simulation results generally fall within one standard error of the physical results.
The response around the anticipated optimum of the broken cross beams model is predicted by our models considerably worse.The fact that the simulations also performed worse in that region was expected, since after cells in the mesh are broken some self-intersection may occur in the simulation.When viewing the best performing straight ITIL sample (broken wb+, Fig. 16(a)) we observe that several failure modes occur throughout the sample (due to manufacturing inaccuracy), but shearing off of part of the cross beams is not one of them.One possible explanation is that for geometry so small as the order of the nozzle size the micro-gaps in between adjacent extruded strands impair the material properties less than in the relatively large test samples with which the material properties of MEX printed parts were determined.The empirically obtained material properties of MEX printed parts may be less accurate at smaller scales.
By comparing the analytical model for the diagonal ITIL variant to the simulations results we observe that at higher   values the analytical model overestimates the strength.The height of the constraint surface for  ta (see Fig. 6(b)) was higher than it was simulated to be, which can be related to the fact that for a higher   the angle  of the beams is lower, which means that the shear component of the stress was higher and that the stress might have been less uniform.
This overestimation of the analytical model is compensated by an underestimation of the simulations, which makes the analytical predictions fit well to the empirical results; see Fig. 14.The underestimation of the simulations can be related to the type of boundary conditions which were applied, but also to an underestimation of the empirically determined material properties of 3D printed parts.Overall the diagonal ITIL model is supported by the data.
The analytical modeling of the dovetail geometries considers only homogeneous tensile stress, but disregards the fact that the interlock is broken by deformation alone.The fact that the dovetail designs do not rely on topological interlocking means that they are more difficult to model.

Comparison of different interlocking structures
All interlocking designs considered can reach roughly 6 to 7 MPa, which means that ±2 MPa of the theoretical upper bound of 8.6 MPa from Eq. ( 1) is used to secure the interlock.The diagonal ITIL design with a cell stress of 7.07 ± 1.0MPa seems to outperform the others, but not significantly.See Fig. 15.The dovetail designs performed considerably worse then the other designs -even when the  max constraint was violated.Moreover, some prints of the dovetail designs showed contamination between the two materials (see Fig. 16(e)), which might be caused by the toolpath discontinuities (see Fig. 11(c)).Some designs of the straight ITIL variant near where the broken model is optimal perform significantly better than the other tested straight ITIL designs; given the high dimensionality of the design space it might be the case that there is a straight ITIL design in between the whole and broken optima which outperforms the diagonal ITIL variant nonetheless.
For designs where a lower  max is allowed the diagonal ITIL variant performs worse.For  max = 1.8 mm the diagonal is simulated to outperform the straight ITIL variant by 5.53 MPa to 4.69 MPa.See Table 2.The performance of the diagonal ITIL model greatly reduces for lower  max , because the angle of the beams is reduced, which causes them to be more susceptible to shear stresses.

Limitations
The biggest obstacle to drawing definitive conclusions about the relative strengths of the interlocking structures is the size of the standard deviation.The variation in results has a multitude of causes, related to print head position inaccuracy, filament diameter fluctuations, temperature oscillations and even the location on the build platform. 2In order to mitigate these factors we have spread out the various geometries over 5 different Ultimaker S5 printers and different locations on the build platform.Although this increases the variability in the results it does increase the reliability.
Because we needed to test a large quantity of samples, we had to limit the size of the samples in terms of interface dimensions.For the 5 × 5 cells of the straight ITIL designs 12 out of 25 cells were boundary cells, so it is challenging to extrapolate these results to larger interfaces.Moreover, the different geometries were tested with different amounts of cells, because the cell geometries were very different.The dovetail designs were modelled without taking friction into account and the size of the interlock (, ) was not optimized for.This makes it hard to draw definitive conclusions comparing the different interlocking geometries.
Finally, the validity of our models turned out to be limited, partly because of homogeneity assumptions and partly because we used empirical material properties of 3D prints, rather than simulating how those properties come about.In order to get more accurate models we could emulate the contact area and polymer entanglement between neighbouring traces, fit that to the empirically obtained tensile properties and use the resulting fitted model to simulate the structures on a toolpath level.

Fig. 1 .
Fig. 1.Applications of dual material products using the ITIL lattice for interlocking Ultimaker transparent PP and Ultimaker green tough PLA material.A gripper making use of the straight ITIL variant.A prosthetic hand using the diagonal ITIL variant.A storage box utilizing the straight ITIL variant.

Fig. 2 .
Fig. 2.Interlocking principles; exploded view with a section cut at the top layer.The dovetail can be disassembled by translation; the cubic lattice causes highly discontinuous extrusion on some layers; ITIL solves these problems.

Fig. 3 .Fig. 4 .
Fig.3.Unit cell of the straight ITIL variant.Force is transferred from material  (green) to material  (transparent cyan) through the contact area on the cross beams.In the partly broken situation there is still interlocking, but the total force is transferred through a narrower area.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.Analytical model for the diagonal ITIL cell.The Z shear constraints are redundant because of the small height of the unit cell.The optimum is determined only by the tensile constraint of PP:  tb .

Fig. 7 .
Fig. 7. Example simulation meshes.The mesh of the diagonal ITIL variant is half a unit cell covering only one of the two diagonal fingers.Because of the rounded corners an axis aligned meshing was impossible, leading to singularities in the mesh.

T
.Kuipers et al.

Fig. 8 .
Fig. 8. 2D slices of the 4D simulation results and fitted RBF hypersurface for the straight ITIL variant.Sampled data points in black.

Fig. 9 .
Fig. 9. Ultimate strength according to the analytical models and the simulation results for the straight variant of the ITIL lattice.The analytical models follow roughly the same shape and same height as the simulation results.The response on three 2D slices of the 3D design space are shown, from left to right at ℎ f = 0.8,   = 2.7 and at   = 2.88.

Fig. 10 .
Fig. 10.Simulation results for diagonal ITIL variant using linear interpolation between the simulation results.

Fig. 11 .
Fig. 11.Gcodes generated with Cura Arachne engine beta.The layers of main fingers are shown and the other layers in a transparent overlay.While the straight and diagonal ITIL lattice produce continuous extrusion beads, the toolpaths for the dovetail designs include small separated segments.

Fig. 12 .
Fig. 12. Experimental setup of straight ITIL variant.Note that for the broken model   is minimal, so va-is not a valid sample.

Fig. 14 .
Fig. 14.Test results compared to the predictions according to the analytical model and the RBF network fitted to the simulation results.The straight ITIL variant samples are labelled relative to the whole and broken optimum, while the rest is labelled by their   value.

Fig. 15 .
Fig. 15.Comparison of the best performing design of each type.The diagonal ITIL design with   = 1.2 mm showed the highest maximum tensile stress.By inspecting the shape of the graph you can determine the dominant failure mode: PLA breaking (sharp drop) or PP yielding (plateau).Dovetail unlocking has either a sharp drop or a more gradual decline.

Fig. 16 .
Fig. 16.Samples after tensile tests of the best performing designs.The broken wb+ samples from the straight ITIL variant exhibit multiple failure modes, indicating that this sample was close to the intersection of several constraint surfaces.

Table 1
Yield properties of single material MEX printed samples.

Table 2
Optimal designs according to the hypersurface fitted to the FEM simulations for the straight ITIL variant and for the diagonal ITIL variant.