 Research
 Open access
 Published:
Measured and simulated mechanical properties of additively manufactured matrixinclusion multimaterials fabricated by material jetting
3D Printing in Medicine volume 10, Article number: 4 (2024)
Abstract
Modern additive manufacturing enables the simultaneous processing of different materials during the printing process. While multimaterial 3D printing allows greater freedom in part design, the prediction of the mixmaterial properties becomes challenging. One type of multimaterials are matrixinclusion composites, where one material contains inclusions of another material. Aim of this study was to develop a method to predict the uniaxial Young’s modulus and Poisson’s ratio of material jetted matrixinclusion composites by a combination of simulations and experimental data.
Fifty samples from commercially available materials in their pure and matrixinclusion mixed forms, with cubic inclusions, have been fabricated using material jetting and mechanically characterized by uniaxial tensile tests. Multiple simulation approaches have been assessed and compared to the measurement results in order to find and validate a method to predict the multimaterials’ properties. Optical coherence tomography and microscopy was used to characterize the size and structure of the multimaterials, compared to the design.
The materials exhibited Young’s moduli in the range of 1.4 GPa to 2.5 GPa. The multimaterial mixtures were never as stiff as the weighted volume average of the primary materials (up to \(33.9\%\) softer for 45% RGD8530DM inclusions in VeroClear matrix). Experimental data could be predicted by finite element simulations by considering a nonideal contact stiffness between matrix and inclusion (\(F_A=2.2\frac{TN}{m^3}\) for RGD8530DM, \(F_A=4.8\frac{TN}{m^3}\) for RGD8430DM), and geometries of the printed inclusions that deviated from the design (rounded edge radii of \(r=220\,\upmu\)m). Not considering this would lead to a difference of the estimation result of up to \(612\,\)MPa (44%), simulating an inclusion volume fraction of 45% RGD8530DM.
Prediction of matrixinclusion composites fabricated by multimaterial jetting printing, is possible, however, requires a priori knowledge or additional measurements to characterize nonideal contact stiffness between the components and effective printed geometries, precluding therefore a simple multimaterial modelling.
Introduction
Additive manufacturing, or 3D printing, systems can do much more than just rapid prototyping. While these systems are often used in industry and the arts to quickly create and improve the design of parts, they are also being used to manufacture parts that would otherwise be difficult or impossible to produce using conventional techniques.
In medicine, additive manufacturing has already reached the point where models can be manufactured directly from the patient image data to train surgeons and to plan complex operations, thus enabling better treatment methods to be developed. Intraoperative surgical tools, as well as implants, have become the new gold standard in medical fields such as orthopedics [1, 2] and dentistry [3, 4].
However, materials ideally suited for medical applications should be as close as possible to natural ones, which come in many different forms and many different properties. Often it is indeed hard to find a singl synthetic material with the same mechanical properties when it comes to replacing or mimicking a human tissue or manufacturing implants and tissue engineering scaffolds [5]. This does not necessarily have to involve materials that are actually used in patient treatment. Most aspiring doctors are trained using body donations, if they even have the opportunity to practice procedures under conditions that approximate reallife scenarios. Since the quality of surgery is significantly influenced by the preparation and training of the surgeon, there is a growing effort to create models for these purposes. Naturally, the closer 3Dprinted materials can mimic real tissue, the better.
Modern additive manufacturing allows to fabricate material composites in great variety, including some that resemble natural materials more closely than others. Mix or multimaterials can be fabricated additively by different means, either by printing two materials close to each other or by mixing the raw materials before depositing them onto the building plate [6].
Among the different multimaterial possibilities, matrixinclusion composites are of special interest for this study. Matrixinclusion composites are materials consisting of one continuous phase, the matrix, containing dispersed inclusions from secondary materials. There are many other possible structures such as discontinuous phase composites or interpenetrating composites with different geometries [7]. As the number of phases increases, so does the complexity and range of possible structures. Human tissues, such as bone, provide examples of how such a composite can be organized into multiple hierarchical levels, each containing phases of different geometries and materials [8,9,10,11]. However, in this study single unconnected inclusions of one phase embedded in a matrix material were considered. Even for this type of composite, many possibilities exist because the shape and distribution of the included phases can be varied.
It is however important to notice that materials can also be mixed during material jet printing. The term digital material has become established for such print products. In this process, two or more primary materials are mixed by the printer to fill one voxel. On the size scale of the print resolution, this results in a homogeneous material [12, 13].
A wide variety of different printing processes to additively manufacture structures, each with their appropriate materials and applications, exists [14, 15]. Multimaterial jetting has been considered here, as this is one key process to fabricate material composites. Hereby photopolymer droplets are dispensed from the print head and subsequently cured with a UV lamp. The resulting mechanical properties of the printed part depend strongly on the used primary materials [16, 17] print orientation [18], and photopolymer curing degree [19,20,21].
Numerous studies have already investigated composites involving micro or nanoparticles, emphasizing the strategic alignment of these particles within the material to achieve specific properties [22,23,24]. In the context of material jetting, however, the inclusions are printed, enabling precise design of their shape and orientation. It should be noted that the inclusions, in this case, are larger than nanoparticles, in a size that is fabricable with the printer resolution (\(>40\,\upmu\)m for \(600\,\)dpi), and also consist of photopolymer similar to the matrix material.
To alter the mechanical properties of the finished part multimaterial printing comes into play, allowing to finetune material mixtures and printed geometries to achieve the desired properties [25, 26]. Prior studies investigated the mechanical properties of multimaterial jetting, by assessing the influence of curing grade and associated hardness of additively manufactured parts [27] or its chemical properties [28]. It has already been shown by Salcedo et al., that test specimens consisting of multiple materials can be simulated by finite element modelling (FEM) [26]. However, only one large inclusion was included in the specimen and strain differences between simulation and experiment were found at the interface of the different materials. This shows the necessity for better modelling of the interfaces, especially in the case of matrixinclusion composite which involves numerous individual inclusions. Other studies investigated estimations of the tensile strength of material jetted multimaterials, where the blending was performed on a voxelwise basis [13, 29, 30]. Various models have been employed to describe mixed material jetted materials, however, none of them deals with a matrixcomposite structure with inclusions ordered in unit cells.
Although there are already approaches to model matrixinclusion composites [31, 32], these methods have not been compared extensively with experimental results from additively manufactured parts. Also, a generalized design approach is missing, to predict and plan multimaterials to achieve specific material properties.
Several studies have reported that the different material phases mix at their interface and form an intermediate layer [26, 29, 30, 33]. The most comprehensive examination was conducted by Zorzetto et al., in which material transitions of material jetted specimens were determined. Using nanoindentation, interface transitions up to \(150\,\upmu\)m in width were identified [33]. Salcedo et al. observed these blurred material transitions both in the scanning electron microscope and the light microscope. However, he did not incorporate them into his simulation model, and his simulation deviates from the experimental results precisely at the material transitions [26]. Bezek et al. utilized the blurred interfaces to blend the two materials voxelwise, aiming to create a smoother transition between the two materials [30]. However, no previous study has investigated how this affects the properties of the overall material and how these effects can be incorporated into a simulation model.
The current study focuses on matrixinclusion composites by using multimaterial jetting. The composite is of a scale that the inclusions can be observed with the bare eye (0.1 to \(1\,\)mm), with the composite becoming nevertheless similar to a homogeneous material in case e.g. of medical models, printed at a scale typically larger than \(10\,\)mm.
In particular, this study examines the fabrication of the samples, their mechanical characterization, and mathematical simulation approaches for predicting composite mechanical properties as a function of the inclusion volume fraction. In addition, the accuracy of the fabricated geometries and the printing process itself are discussed in terms of their implications for improving the predictive ability of simulations.
This article is structured such that it begins by elucidating the methods used for 3D printing the test specimens. Followed by a description of the geometric characterization and the mechanical testing procedure. Subsequently, the two simulation methods are explained. In the next section, the results about to the aforementioned steps are presented, before concluding with a discussion of these findings.
Materials and methods
Samples were produced with an Connex3 Objet500 polyjet printer (Stratasys Ltd., Minnesota, USA), from the raw materials VeroClear, VeroPureWhite, TangoPlus and TangoBlackPlus, all of them are acrylicbased photopolymers. The printer was operated with the software Objet Studio Ver. 9 (Stratasys Ltd., Minnesota, USA). The material jet printer has a lateral resolution of \(600\,\)dpi and used a layer thickness of \(30\,\upmu\)m, therefore the voxels have a dimension of \(42.3\times 42.3\times 30\,\upmu\)m.
After the printing, the samples have been cleaned and dried. All samples have been processed with the socalled “matte” setting, covering the whole sample in support material and therefore producing a uniform surface on the top and bottom side of the sample. Alternatively, there is the “glossy” setting, where no support material is applied, resulting in different surface structures on the sample’s upper side and the bottom side facing the build platform. However, this machine setting does not affect the interfaces between matrix and inclusions.
Sample sets of 5 specimens each are considered, according to ISO 527 [34]. By fabricating VeroPureWhite together with TangoBlackPlus, the digital material RGD8530DM can be created. Multimaterial samples with 3 different inclusion volume fractions \(f_{inc}=\)10%, 30% and 45% were designed using VeroClear as the matrix material and RGD8530DM for the inclusions. From preliminary prints and literature is known, that the material phases combine with a certain layer thickness at the interface. Therefore samples with higher volume fractions have not been printed to avoid unintentional connection of the inclusions, since the aim of the study was to investigate the matrixinclusion composites. Also, cubical inclusions were chosen, since other inclusion geometries would connect to each other at lower volume fractions. For example, spheres would have a diameter of \(1\,\)mm at a volume fraction of \(f=52.4\%\) in a \(1\,\)mm\(^3\) cubical unit cell.
The 3 sets of multimaterials have been printed together with the 2 sets of pure primary materials in a single print job. All samples were printed with the same orientation, the long axis of the specimen in the direction of the print head movement.
A second print job was done using the same setup as stated before, but with RGD8430DM as the inclusion material, which is fabricated using VeroPureWhite and TangoPlus. Altogether 10 sets with a total of 50 samples were prepared and tested.
These materials were selected based on their established miscibility and adhesive characteristics at the inclusionmatrix interface. For higher contrast in the OCT and visibility of the inclusions in the microscope, a transparent matrix and opaque inclusion materials were chosen. All samples were stored and dried for three days between printing and mechanical testing.
Specimen geometry
The specimen geometry used is the “Type A” multipurpose test specimen from ISO 3167 [35], as required by ISO 527 [34], often referred to as dogbone. All specimens contain inclusions in a regular pattern, arranged in cubic unit cells, where the inclusion is in the center of the unit cell. Each unit cell has a side length of 1 mm, consisting of matrix material on the outside and an inclusion in its center. The placement in the tensile specimen was done in a way, that in the cross section of the gauging area fitted \(10\times 4\) complete unit cells. The dimensions of the inclusions have been chosen in a way, that the inclusions account for a certain proportion of the total volume, for example, a cubic inclusion with an edge length of of \(464\,\upmu\)m would have a volume fraction of 10% in a \(1\,\textrm{mm}^3\) unit cell.
Measurements
The measurement of the mechanical properties was performed according to ISO 527 [34]. The specimens were fixed in the test stand Messphysik BETA 102,5 (Messphysik GmbH, Fürstenfeld, Austria) and then pulled apart with a feed rate of \(0.75\,\frac{\textrm{mm}}{\textrm{min}}\), which was increased to \(3\,\frac{\textrm{mm}}{\textrm{min}}\) after 1.5 mm of displacement. The pretension was 1N and the length and width of the specimen were observed with an OS65D video extensometer (Mintron Enterprise Co., Ltd., Taipei, Taiwan). The values were captured at \(50\,\)Hz. The Young’s modulus was evaluated at strains from \(\varepsilon _1 = 0.05\%\) to \(\varepsilon _2 = 0.25\%\) in the linear range of the material from the measured force. The start and end values of the measurement were determined by linear regression in the measuring range. For lateral strain, a 20frame moving average was applied before calculating Poisson’s ratio to reduce scatter.
To characterize the dimensions of inclusions in the printed specimen, three optical methods have been employed: optical microscopy and profilometry as well as optical coherence tomography. Microscope images were taken with the Levenhuk DTX 90 (Levenhuk, Inc., Florida, USA) and VHX7000 (Keyence International, Mechelen, Belgium) digital microscopes. Profilometer measurements of the sample surface were made with the optical profilometer Keyence VR5000 (Keyence International, Mechelen, Belgium).
The samples have been inspected with a custommade sweptsource optical coherence tomography (OCT) system operating at 1300 nm with a bandwidth of 30 nm [36]. This wavelength allows to easily image through the whole phantom thickness. A field of view (FoV) of \(10\times 10\,\)mm was acquired, the Bscan direction being roughly perpendicular to the printing direction of the phantom. After standard OCT postprocessing, the 3D images were rotated to have the inclusion pattern coalign with the imaging axis. An average intensity projection over a depth of 812\(\,\upmu\)m (see Fig. 7), covering the first inclusion layer, was used to estimate the edge radius of the inclusions.
For a closer look at the microstructure, a second, microscopic spectral domain OCT was used, working at \(800\,\)nm with a bandwidth of 200 nm [37].
Simulation
Two approaches for the simulation of the mechanical behavior of multimaterial properties were used: numerical finite element analysis (finite element method, FEM) was performed with ANSYS (ANSYS, Inc, Canonsburg, USA) and analytical continuum mechanical simulations, derived from the MoriTanaka scheme [38], performed with MATLAB (The MathWorks, Inc., Natrick, USA).
Finite element analysis
Although additively manufactured parts have anisotropic mechanical properties [18, 39, 40], the FEM simulations have been performed with an isotropic material model as an approximation for the uniaxial tensile tests. The primary material properties (Young’s modulus and Poisson’s ratio) have been obtained from the measurements. The volume fraction of the inclusion was varied in \(5\%\) steps. Input for the design was the side length of the inclusion, while the edge radius was considered to achieve the correct volume fractions. To model the interfaces between inclusion and matrix, a bonded contact with an absolute value for normal stiffness was provided.
Mimicking the test setup, the FEM analysis was performed with \(10\times 4\) unit cells of one millimeter cell length, which equates to one layer of unit cells within the tensile test specimen, as shown in Fig. 1. There are no interfaces between the individual unit cells since all matrix phases have been merged to one single geometry. The geometries were created with the ANSYS tool “DesignModeler”, the definition of boundary conditions, meshing and simulation were done with the ANSYS tool “Mechanical”.
A strain of \(\varepsilon =0.25\%\) was applied transverse to the unit cell plane, similarly to the range applied during experimental testing, while the other side was fixated in the direction along the strain. The average stress in the direction of the elongation has been recorded to calculate the macroscopic Young’s modulus from the simulated uniaxial stress and the preset strain. The Poisson’s ratio of the material mixtures was evaluated by measuring the average displacement of the geometrical side planes compared to the preset strain. The calculation of the material parameters from this simulation setup was validated with pure materials.
Adaptive meshing size was used with a maximum element size of \(300\,\upmu\)m, using the finest resolution from the standard meshing tool, using about \(2.1\cdot 10^6\) to \(5.9\cdot 10^6\) mesh nodes and \(1.4\cdot 10^6\) to \(3.7\cdot 10^6\) mesh elements, depending on the inclusion geometry. The triangular meshes were checked for convergence. A higher number of mesh nodes was required for inclusions with very high volume fraction (\(f_{inc}\rightarrow 85\%\)) due to the larger material interface and for very small inclusions (\(f_{inc}\rightarrow 5\%\)) because of the potential for more movement in the matrix. The contact between matrix and inclusions was defined as bonded, but with different normal stiffnesses \(F_A\). The average stress in the direction of the elongation has been recorded to calculate the Young’s modulus. Three values for this stiffness were evaluated as well as idealized bonding (\(F_A\rightarrow \inf\)).
Degeneration of the inclusion geometry affects cubic inclusions (e.g. rounding of edges) during the printing process. Therefore this geometric difference has also been simulated, by adding fillets to the inclusion edges with different radii.
Parameter optimization was performed to find the ideal values for the edge radius and contact stiffness. Simulations were performed with different edge radii (step size \(10\,\upmu\)m) and contact stiffness (step size \(100\,\frac{\textrm{GN}}{\textrm{m}^3}\)). The parameter set with the lowest accumulated deviation from the mean of the three measured sample sets was selected.
Inclusions with high volume fractions (\(f_{inc}>85\%\)), although not measured experimentally, were also simulated. Here the inclusions are interconnected and not separated within the matrix. The required volume fraction to connect the inclusions depends on their edge radius. Simulation steps beyond this value were chosen with a step width of \(2.5\%\) of the inclusion volume.
Multimaterial homogenization
Multimaterial homogenization was performed to evaluate whether a model of ideal inclusion geometries (cubes and spheres), which does not consider inclusionmatrix interface properties, is still able to reproduce experimental findings. To parameterize this model the elastic compliance \(\textbf{C}\), evaluated from the measurement values, was represented as a secondorder tensor in KelvinMandel notation. Secondorder tensor representation was used for the Matlab implementation. The simulation was performed according to the MoriTanaka scheme, where inclusions are assumed to be in a homogeneous matrix material [38, 41, 42]. The volume increment of the simulation was set at \(1\%\). The simulation models mixing scheme was adapted from Fritsch and Hellmich [8] and the Eshelby tensors from Yozo Mikata [43, 44]. For more details on multimaterial homogenization and the corresponding Hill and Eshelby tensors in the implemented form, please refer to Appendix A.
Results
No misprints occurred during production. In the transparent matrix, opaque inclusions are visible with the bare eye, as shown in Fig. 2, and the different inclusion volume fractions can be distinguished. In the direction, normal to the 3D printing plane, the transparent matrix allows to look through the material and gives a good view on the inclusions, however, traces from the printing are visible under the microscope. The side surfaces of the samples are much rougher than the top and bottom surfaces, and although inclusions are still visible, the otherwise transparent materials appear almost opaque.
Tensile tests
The measured uniaxial Young’s modulus for the samples of the pure primary material RGD8530DM was \(1519.1\pm 39.9\,\)MPa, and for RGD8430DM it was \(1773.8\pm 29.9\,\)MPa. The two sample sets of pure VeroClear, originating from the two print jobs with different mixmaterials, yielded Young’s moduli of \(2431.0\pm 29.6\,\)MPa (printed together with RGD8430DM samples) and \(2444.6\pm 36.1\,\)MPa (printed together with RGD8530DM samples), with an insignificant deviation (twosample ttest pvalue for nonequal means \(p = 0.5328\)).
All tested multimaterial samples have a Young’s modulus below the volumeaveraged mean of the primary materials. The multimaterials show a lower Young’s modulus than pure VeroClear samples, which further decreases with a higher content of the softer inclusion material. The results with RGD8530DM and RGD8430DM inclusions are plotted in Fig. 3. The detailed list of measured values is given in Appendix B.
Representative curves showing the technical stress over strain in samples from the print with RGD8530DM are shown in Fig. 3c. The multimaterials rupture at earlier elongations than the samples made from pure materials. While the multimaterials ruptured at strains ranging from \(2.3\%\) to \(4.5\%\), the pure materials ruptured at strains ranging from \(6.2\%\) to \(22.7\%\), with the highest strains occurring in the pure RGD8430DM samples.
The samples with inclusions tore at the interface between inclusion and matrix material, indicating this cross section as the weakest part of the multimaterial. The height profile of a typical rupture area is shown in Fig. 4.
OCT images
Since the used materials have a similar radiodensity, only little contrast can be achieved by computer tomography, but the inclusions with different colors are distinguishable by OCT imaging. For further image processing, planes of the materials containing inclusions have been investigated, and the selected region of interest (ROI) is the first layer of inclusions. Using OCT measurements with pure materials, a refractive index of \(n=1.55\) for VeroClear and \(n=1.5175\) for RGD8530DM was found.
Shown in Fig. 5 is a projection of the volume scan, where the cubical inclusions can be seen in the matrix material. However, these cubes do not have a perfectly rectangular projection, but the rounded edges can be seen. Rounded edges with radii up to \(r\!=\!400\,\upmu\)m were found, while most edges have a radius from \(r=200\) to \(250\,\upmu\)m.
From the image data, the planes of inclusions have been stacked to a single image and have been segmented to distinguish the inclusions and matrix more easily. The designed inclusion edge length for \(1\!\times \!1\!\times \!1\,\)mm unit cells with a volume fraction of \(f_{inc}=30\%\) is \(s_{inc}=669\,\upmu\)m. Measurements showed that the cubes were smaller transversal to the printing direction, up to \(\Delta l_p\!=\!15\,\upmu\)m. But the inclusions were larger, in the printing direction, up to \(\Delta l_t\!=\!60\,\upmu\)m. For comparison, the material jetting 3D printer has a resolution of \(600\,\)dpi, meaning that each voxel has a side length of \(42.3\,\upmu\)m.
In Fig. 6 the inclusions are shown and the individual layers can be distinguished. It is important to notice that OCT images do not represent precise geometric distances, but instead optical distances that depend on the refractive indices of the materials and their transitions. However, these scans convey the information that between the layers there is a sharp transition between inclusions and matrix, but within the layers, the materials are more fringed, which can lead to different interface properties.
Also shown in Fig. 6 are the same inclusion, but with an orthogonal profile. A geometric deformation is visible, with a wavy surface on each cube. The waviness increases as more layers are formed during the printer’s operation. Not visible in this image is that the waviness of the layers decreases as the inclusions stop. After a few layers consisting purely of matrix material, the surface is planar again and the next inclusions start from flat again.
Simulations
Figure 7 shows the fit of the simulation of the matrixinclusion composite to the measurement results. For RGD8530DM inclusions in a VeroClear matrix, simulation results obtained using a contact stiffness of \(2.2\frac{\textrm{TN}}{\textrm{m}^3}\) and a radius of \(r=220\,\upmu\)m that matched the OCT images, are within the interquartile range of the measurements. Accordingly, the inclusions will be in contact with each other at a volume fraction of \(f_{inc,th}=88.96\%\) for \(r=220\,\upmu\)m. FEM simulations with ideal matrixinclusion bonds are also shown and closely match the results from the multimaterial homogenization for perfectly cubic and spherical inclusions. The expected Young’s moduli from multimaterial homogenization are lower, but close, to the volume averaged mean of the pure materials (black dashed reference line), while the FEM simulation with ideally bonded contacts gives both slightly lower and higher values to the volume averaged mean.
Stiffer contact conditions between inclusion and matrix almost always mean that the entire material is also stiffer, with a translation of the overall stiffness at each volume inclusion ratio as shown in Fig. 8a. Also, a slight change in the inclination of the slope is observable. However, a change in the edge radii of the inclusions means that the behavior of the material mixtures differs with the curvature of the stiffness vs. inclusion fraction changing with changing radii Fig. 8b.
For multimaterials with RGD8430DM inclusions, the same procedure was performed as well. According to the parameter optimization study, simulations also fit with an edge radius \(r=220\,\upmu\)m, however, a higher contact stiffness of \(4.8\frac{\textrm{TN}}{\textrm{m}^3}\) was necessary to ensure a good fit.
Poisson’s ratio
The measurement of Poisson’s ratio using the video extensometer provided measurements with high variability with standard deviations up to 0.08 (valid values for Poisson’s ratio \(0\le \nu \le 0.5\)). Therefore the median values of the pure materials were used as an input for the model. The values can be found in Appendix B Table 1.
The effects of different Poisson’s ratios were investigated by simulation. The results from the multimaterial homogenization show an arced behavior from one of the assumed Poisson’s ratios, known from this measurement or literature, for VeroClear to the Poisson’s ratio of RGD8430DM \(\nu _{RGD8430}=0.415\) and RGD8530DM \(\nu _{RGD8530}=0.418\). Depicted in Fig. 9 are the results for multimaterial homogenization with cubical inclusions in full lines and with spherical inclusions in dashed lines.
The Poisson’s ratios for the simulated multimaterials with ideally bonded interfaces and rounded cubical inclusions are in between the values of simulations with perfectly cubic or spherical inclusions. FEM simulations with either perfect cubic or spherical inclusions gave nearly the same results as the associated multimaterial homogenization (differences of the Poisson’s ratio \(\ll 0.01\)). Contrary to this, the simulations with springloaded interfaces show a decreasing Poisson’s ratio for higher inclusion fractions, away from the Poisson’s ratio of the inclusion material.
Discussion
Results have shown that the material jetted multimaterials have Young’s moduli lower than the weighted average of the primary materials, in some cases even lower than that of the softer primary material. The relatively small standard deviations speak for the reliability of the measurements. Simulation results showed that welltested methods to predict the properties of multimaterials, like multimaterial homogenization, cannot be applied easily, since the process of 3D printing introduces novel effects that influence mechanical material behavior.
While the results of the multimaterial homogenization agree with the FEM simulation with ideally bonded interfaces, a closer look at the samples (Fig. 2) reveals a zone where matrix and inclusion material is mixed. This mixture between the two material phases is most likely the cause of the difference between the two methods, and most likely the cause of the discrepancies between measurements and simulations under the assumption of ideal interface contacts. The fringed character of the interfaces can indeed be observed in the OCT images and confirms previous observations [33, 45]. The interface within a layer appears blurry, most likely because the jetted material droplets have a short time to blend before they are cured by UV light. This is also assumed the reason, why the cubes have radii on the edges, larger than the printer resolution.
Simulations with a finite contact stiffness, representing a springloaded interface, yield results with lower Young’s moduli. These simulations match better the measurement results, and choosing the right contact stiffness allows the prediction of the multimaterials properties. Especially when also the inclusion deformation (rounded edges) has been taken into account, simulations with a good fit to the measurement results could be established. The simulations show that the inclusion interface is essential to improve multimaterial predictions, and what the effect of varying inclusion geometries can be.
Although the simulation with a specific contact stiffness yields good results for Young’s modulus, it shows an unexpected behavior for the Poisson’s ratio, since this value is decreasing for higher inclusion contents, instead of rising to the Poisson’s ratio of the inclusion primary material. This is most likely due to the implementation of normal stiffness in the simulation software, as pulling on a spring, which represents the interfaces, is separating the surfaces of the inclusion and matrix. This causes a displacement that is not accounted for in the width change of the simulated geometry and therefore distorts the Poisson’s ratio. Possibly this model can be improved by implementing different stiffnesses, depending on the applied load and direction. For example finite stiffness values for tension, but infinity stiffness values (\(F_a \rightarrow \inf\)) for compression.
However, both the multimaterial homogenization, as well as FEM simulation with ideally bonded interfaces yield Poisson’s ratios whose trends coincide with the measured values. But changing the Poisson’s ratio in the simulation has only a minor impact on the Young’s modulus, which is in most applications the more important property to consider.
Since the literature provides many different values for the Poisson’s ratio of the primary material VeroClear [7, 46, 47], additional simulations, not shown in this article, were performed, where the Poisson’s ratio assumed for VeroClear was varied between the following values \(\nu =[0.33;\, 0.35;\, 0.38;\, 0.41]\) in addition to the measured values of the pure materials. The Poisson’s ratio of VeroClear has a different impact on the simulation with ideally bonded interfaces (for rounded cubical RGD8530DM inclusions: difference between \(\nu _{VC}=0.33\) and \(\nu _{VC}=0.41\): \(13.9\,\)MPa at \(f_{inc}=10\%\); \(13.0\,\)MPa at \(f_{inc}=50\%\); \(2.17\,\)MPa at \(f_{inc}=85\%\)), than on the simulation with a surface stiffness of \(F_A = 2.2\,\frac{\textrm{TN}}{\textrm{m}^3}\) (difference between \(\nu _{VC}=0.312\) and \(\nu _{VC}=0.41\): \(4.34\,\)MPa at \(f_{inc}=10\%\); \(18.9\,\)MPa at \(f_{inc}=50\%\); \(21.2\,\)MPa at \(f_{inc}=85\%\)).
While the difference of the multimaterials Young’s moduli, depending on the matrix materials Poisson’s ratio, decreases at higher volume fractions for simulations with ideally bonded contacts, the simulations with springloaded interface conditions show the exact opposite behavior. However, for materials in the range of about \(1\,\)GPa, the change of a few MPa is insignificant.
Knowing the properties of one measured multimaterial, with a specific volume fraction of the inclusions, is sufficient to interpolate the contact conditions for a relatively simple model as outlined in this study, and therefore similar multimaterials with differing volume fractions can be predicted. However, we are still lacking a material model to predict mixed materials just from the properties of the primary materials, since we do not know how the interface conditions will turn out for new materials.
This can be seen particularly well when comparing the multimaterials with RGD8530DM and RGD8430DM inclusions. Here, the inclusion geometry depends mainly on the material jetting process, so the geometry, especially the edge radius, is similar. However, these different materials bond differently to the matrix material, therefore the contact stiffness varies.
Still, there are many further parameters that can be changed in FEM simulations, for instance, the tangential stiffness of the interface, since only variations in the normal stiffness were simulated in this study. Even more, when assuming an intermediate material between the matrix and the inclusion, there are suddenly many more parameters to be defined, like the position, geometry and contact conditions of the intermediate material, and of course the material properties itself. While it is easy to measure the properties of the primary materials, unfortunately, the contact stiffness or properties of an intermediate layer cannot be characterized with the same standardized tensile test, which in turn means, that these properties have to be estimated from indirect measurements. To emulate different interface properties, which are common to polyjetted materials, and especially pronounced with cubic inclusions, the boundary conditions for the contacts in layer (four lateral sides of the inclusions) can be assumed different to the conditions to the next layer (top and bottom sides).
In the field of multimaterial homogenization, there are also some approaches to model inclusions consisting of multiple layers [31, 32, 48, 49]. A possible improvement could be to represent the contact properties as an encapsulation of the inclusions.
Returning to the geometric deviations, OCT observed not only edge rounding but also nonflat, but wavy layers. Although it was beyond the scope of this study to account for this variation, it would be useful to investigate this effect in more detail. In particular, it should be determined whether the waviness of the layers in the printing direction occurs in all printing models, or whether it is only a peculiarity of this particular printing device. Depending on this, it may also be possible to determine under which conditions and geometries waviness occurs. If it turns out to be a common phenomenon, eventually caused by different shrinkage rates of the materials during photopolymerization [18], it will be necessary to investigate how such a geometry deviation can be predicted for multimaterials and considered in the material models.
Further observations from this and preliminary studies, which need further investigation, are the shelf time of the material cartridges. If the material even slightly changes its properties in the cartridge, although it is still usable for 3D printing, the resulting elasticity might change. This would have an impact on the repeatability and predictability of the material properties that have to be evaluated.
Also storing the materials after the print has not been investigated, as the samples were tested shortly after the production. However, it is quite conceivable that in applications where the elasticity of the material plays an important role, it should remain unchanged over a long period of time. Material aging, in particular the effect of composite materials and their contact points, still needs to be investigated.
It is reasonable to assume that this simulation concept applies to all polymers produced through material jetting. Since there are various highly elastic, rubberlike printing materials available, this expands the range of printable multimaterials further. Naturally, the parameters for contact stiffness would differ for other primary materials. Additionally, this FEM approach should work for other composites, regardless of how the individual components are arranged within the volume or whether more than two individual materials were used, since the geometry does not change the concept of nonideal contact stiffness. A different arrangement of material phases will alter the overall behavior of the 3Dprinted component. Especially concerning anisotropy, the freedom to choose inclusion geometries provides ample room for material design. This, in turn, could be highly valuable for creating as realistic as possible models for educational purposes in the medical field or for surgical planning.
An alternative arrangement of inclusions in a matrixinclusion composite could also offer advantages in terms of tensile strength. This is because the experiments have shown that the interfaces between inclusions and the matrix are the weakest points of the material. By placing the inclusions in a way that eliminates the presence of a fracture plane, it becomes possible to achieve different fracture behaviors.
Conclusion
In this study, material jetted matrixinclusion composites, where the individual phases can be distinguished with the bare eye, were investigated for their mechanical and geometrical properties. This is the first time that this type of material jetted compound has been studied in detail, and conclusions drawn from this study can be applied to other geometries. Comprehensive material tests and OCT examinations were able to reveal some peculiarities of the polyjet printer, including deviations from the planned design. It has been observed that multimaterials depend not only on the starting materials, but also very much on the geometries and contact properties between the individual phases. This work confirms the findings that transition layers form between material phases, which have also been reported in several other articles. However, this is the first time that the effects of these deviations on material properties have been studied and incorporated into a simulation model.
For this kind of multimaterial, the Young’s modulus is not trivial to predict. Tuning the interface conditions between inclusion and matrix in the FEM simulation allows to estimate the properties of the resulting additively manufactured material. But these FEM simulations have a great variety of parameters to set, which already require knowledge of the multimaterials behavior. Since, to our best knowledge, no precise models exist, to predict the interface properties of 3D printed multimaterials, a minimum of empiric knowledge still is required as input for this model to achieve the desired material properties.
Availability of data and materials
Data will be made available on request.
References
Wong KC. 3Dprinted patientspecific applications in orthopedics. Orthop Res Rev. 2016;8:57–66. https://doi.org/10.2147/ORR.S99614.
Liu K, Zhang Q, Li X, Zhao C, Quan X, Zhao R, et al. Preliminary application of a multilevel 3D printing drill guide template for pedicle screw placement in severe and rigid scoliosis. Eur Spine J. 2017;26(6):1684–9. https://doi.org/10.1007/s0058601649261.
Bai L, Zheng L, Ji P, Wan H, Zhou N, Liu R, et al. Additively Manufactured Latticelike Subperiosteal Implants for Rehabilitation of the Severely Atrophic Ridge. ACS Biomater Sci Eng. 2022;8(2):912–20. https://doi.org/10.1021/acsbiomaterials.1c00962.
Pillai S, Upadhyay A, Khayambashi P, Farooq I, Sabri H, Tarar M, et al. Dental 3DPrinting: Transferring Art from the Laboratories to the Clinics. Polymers. 2021;13(1):157. https://doi.org/10.3390/polym13010157.
Melchels FPW, Feijen J, Grijpma DW. A review on stereolithography and its applications in biomedical engineering. Biomaterials. 2010;31(24):6121–30. https://doi.org/10.1016/j.biomaterials.2010.04.050.
Han D, Lee H. Recent advances in multimaterial additive manufacturing: methods and applications. Curr Opin Chem Eng. 2020;28:158–66. https://doi.org/10.1016/j.coche.2020.03.004.
Su FY, Sabet FA, Tang K, Garner S, Pang S, Tolley MT, et al. Scale and size effects on the mechanical properties of bioinspired 3D printed twophase composites. J Mater Res Technol. 2020;9(6):14944–60. https://doi.org/10.1016/j.jmrt.2020.10.052.
Fritsch A, Hellmich C. ‘Universal’ microstructural patterns in cortical and trabecular, extracellular and extravascular bone materials: Micromechanicsbased prediction of anisotropic elasticity. J Theor Biol. 2017;244(4):597–620. https://doi.org/10.1016/j.jtbi.2006.09.013.
Christian H, FranzJosef U. Micromechanical Model for Ultrastructural Stiffness of Mineralized Tissues. J Eng Mech. 2002;128(8):898–908. https://doi.org/10.1061/(ASCE)07339399(2002)128:8(898).
Hellmich C, Barthélémy JF, Dormieux L. Mineralcollagen interactions in elasticity of bone ultrastructure  a continuum micromechanics approach. Eur J Mech A/Solids. 2004;23(5):783–810. https://doi.org/10.1016/j.euromechsol.2004.05.004.
Hellmich C, Ulm FJ, Dormieux L. Can the diverse elastic properties of trabecular and cortical bone be attributed to only a few tissueindependent phase properties and their interactions? Biomech Model Mechanobiol. 2004;2(4):219–38. https://doi.org/10.1007/s1023700400400.
Ituarte IF, Boddeti N, Hassani V, Dunn ML, Rosen DW. Design and additive manufacture of functionally graded structures based on digital materials. Additive Manuf. 2019;30:100839. https://doi.org/10.1016/j.addma.2019.100839.
Yu KM, Tang YM, Chan LC. Property estimate for inkjet based direct digital manufacturing. Comput Aided Des Appl. 2018;15(2):188–92. https://doi.org/10.1080/16864360.2017.1375668.
Jandyal A, Chaturvedi I, Wazir I, Raina A, Ul Haq MI. 3D printing  A review of processes, materials and applications in industry 4.0. Sustain Oper Comput. 2022;3:33–42. https://doi.org/10.1016/j.susoc.2021.09.004.
Tee YL, Peng C, Pille P, Leary M, Tran JP. PolyJet 3D Printing of Composite Materials: Experimental and Modelling Approach. JOM J Miner Met Mater Soc. 2020;72:13. https://doi.org/10.1007/s1183702004014w.
Slapnik J, Pulko I. Tailoring properties of photopolymers for additive manufacturing with mixture design. Progress Additive Manuf. 2021;6(1):83–91. https://doi.org/10.1007/s40964020001471.
Ligon S, Schwentenwein M, Gorsche C, Stampfl J, Liska R. Toughening of photocurable polymer networks: A review. Polym Chem. 2015;7. https://doi.org/10.1039/C5PY01631B.
Königshofer M, Stoiber M, Unger E, Grasl C, Moscato F. Mechanical and Dimensional Investigation of Additive Manufactured Multimaterial Parts. Front Phys. 2021;9. https://doi.org/10.3389/fphy.2021.635736.
Bennett J. Measuring UV curing parameters of commercial photopolymers used in additive manufacturing. Additive Manuf. 2017;18:203–12. https://doi.org/10.1016/j.addma.2017.10.009.
Hofstetter C, Orman S, Baudis S, Stampfl J. Combining cure depth and cure degree, a new way to fully characterize novel photopolymers. Additive Manuf. 2018;24:166–72. https://doi.org/10.1016/j.addma.2018.09.025.
Fuh JYH, Lu L, Tan CC, Shen ZX, Chew S. Curing characteristics of acrylic photopolymer used in stereolithography process. Rapid Prototyp J. 1999;5(1):27–34. https://doi.org/10.1108/13552549910251855.
Niendorf K, Raeymaekers B. Additive Manufacturing of Polymer Matrix Composite Materials with Aligned or Organized Filler Material: A Review. Adv Eng Mater. 2021;23(4):2001002. https://doi.org/10.1002/adem.202001002.
Xu W, Jambhulkar S, Ravichandran D, Zhu Y, Kakarla M, Nian Q, et al. 3D PrintingEnabled Nanoparticle Alignment: A Review of Mechanisms and Applications. Small. 2021;17(45):2100817. https://doi.org/10.1002/smll.202100817.
Xu W, Jambhulkar S, Zhu Y, Ravichandran D, Kakarla M, Vernon B, et al. 3D printing for polymer/particlebased processing: A review. Compos Part B Eng. 2021;223:109102. https://doi.org/10.1016/j.compositesb.2021.109102.
Doubrovski EL, Tsai EY, Dikovsky D, Geraedts JMP, Herr H, Oxman N. Voxelbased fabrication through material property mapping: A design method for bitmap printing. ComputAided Des. 2015;60:3–13. https://doi.org/10.1016/j.cad.2014.05.010.
Salcedo E, Baek D, Berndt A, Ryu JE. Simulation and validation of three dimension functionally graded materials by material jetting. Additive Manuf. 2018;22:351–9. https://doi.org/10.1016/j.addma.2018.05.027.
Yang Y, Li L, Zhao J. Mechanical property modeling of photosensitive liquid resin in stereolithography additive manufacturing: Bridging degree of cure with tensile strength and hardness. Mater Des. 2019;162:418–28. https://doi.org/10.1016/j.matdes.2018.12.009.
Safka J, Ackermann M, Martis D. Chemical resistance of materials used in additive manufacturing. MM Sci J. 2018. https://doi.org/10.17973/MMSJ.2016_12_2016185.
Kaweesa D, Bobbio L, Beese AM, Meisel NA. Prediction and validation of composite mechanical properties resulting from voxelbased microstructural design in material jetting. Rapid Prototyp J. 2023;29(3):488–503. https://doi.org/10.1108/RPJ1220210343.
Bezek LB, Chatham CA, Dillard DA, Williams CB. Mechanical properties of tissuemimicking composites formed by material jetting additive manufacturing. J Mech Behav Biomed Mater. 2022;125:104938. https://doi.org/10.1016/j.jmbbm.2021.104938.
Vieville P, Bonnet A, Lipinski P. Modelling Effective Properties of Composite Materials Using Inclusion Concept. General Considerations. Arch Mech. 2006;58:201–39.
Herve E, Zaoui A. nLayered inclusionbased micromechanical modelling. Int J Eng Sci. 1993;31(1):1–10. https://doi.org/10.1016/00207225(93)900594.
Zorzetto L, Andena L, BriaticoVangosa F, De Noni L, Thomassin JM, Jérôme C, et al. Properties and role of interfaces in multimaterial 3D printed composites. Sci Rep. 2020;10(1):22285. https://doi.org/10.1038/s41598020792300.
Standard DIN EN ISO 5271:2019. Plastics  Determination of tensile properties  Part 1: General principles. 2019. https://doi.org/10.31030/3059426.
Standard DIN EN ISO 3167:2003. Plastics  Multipurpose test specimens. ICS 83.080.01 (International Classification for Standards). 2003.
Krainz L, Veselka L, Mindrinos L, Drexler W, Elbau P. A Quantitative Model for Optical Coherence Tomography. Sensors. 2021;21(20). https://doi.org/10.3390/s21206864.
Giardina G, Micko A, Bovenkamp D, Krause A, Placzek F, Papp L, et al. MorphoMolecular Metabolic Analysis and Classification of Human Pituitary Gland and Adenoma Biopsies Based on Multimodal Optical Imaging. Cancers. 2021;13(13). https://doi.org/10.3390/cancers13133234.
Mori T, Tanaka K. Average stress in matrix and average elastic energy of materials with misfitting inclusions. Acta Metall. 1973;21(5):571–4. https://doi.org/10.1016/00016160(73)900643.
Lumpe TS, Mueller J, Shea K. Tensile properties of multimaterial interfaces in 3D printed parts. Mater Des. 2019;162:1–9. https://doi.org/10.1016/j.matdes.2018.11.024.
Torrado A, Roberson D. Failure Analysis and Anisotropy Evaluation of 3DPrinted Tensile Test Specimens of Different Geometries and Print Raster Patterns. J Fail Anal Prev. 2016;16. https://doi.org/10.1007/s1166801600674.
Zaoui A. Structural Morphology and Constitutive Behaviour of Microheterogeneous Materials. In: Suquet P, editor. Continuum Micromechanics. Vienna: Springer Vienna; 1997. p. 291–347. https://doi.org/10.1007/9783709126622_6.
Zaoui A. Continuum Micromechanics: Survey. J Eng Mech. 2002;128(8):808–16. https://doi.org/10.1061/(ASCE)07339399(2002)128:8(808).
Mikata Y. Determination of piezoelectric Eshelby tensor in transversely isotropic piezoelectric solids. Int J Eng Sci. 2000;38:605–41. https://doi.org/10.1016/S00207225(99)000506.
Mikata Y. Explicit determination of piezoelectric Eshelby tensors for a spheroidal inclusion. Int J Solids Struct. 2001;38(40):7045–63. https://doi.org/10.1016/S00207683(00)004194.
De Noni L, Zorzetto L, BriaticoVangosa F, Rink M, Ruffoni D, Andena L. Modelling the interphase of 3D printed photocured polymers. Compos Part B Eng. 2022;234:109737. https://doi.org/10.1016/j.compositesb.2022.109737.
Javid F, SmithRoberge E, Innes MC, Shanian A, Weaver JC, Bertoldi K. Dimpled elastic sheets: a new class of nonporous negative Poisson’s ratio materials. Sci Rep. 2015;5(1):18373. https://doi.org/10.1038/srep18373.
Donahue CM, Anzel PWJ, Bonanomi L, Keller TA, Daraio C. Experimental realization of a nonlinear acoustic lens with a tunable focus. Appl Phys Lett. 2014;104(1):014103. https://doi.org/10.1063/1.4857635.
Scheiner S, Komlev VS, Gurin AN, Hellmich C. Multiscale Mathematical Modeling in Dental Tissue Engineering: Toward ComputerAided Design of a Regenerative System Based on Hydroxyapatite Granules, Focussing on Early and MidTerm Stiffness Recovery. Front Physiol. 2016;7. https://doi.org/10.3389/fphys.2016.00383.
Bertrand E, Hellmich C. Multiscale Elasticity of Tissue Engineering Scaffolds with TissueEngineered Bone: A Continuum Micromechanics Approach. J Eng Mech. 2009;135(5):395–412. https://doi.org/10.1061/(ASCE)07339399(2009)135:5(395).
Acknowledgements
The authors would like to thank Gunpreet Oberoi for her valuable help in the preparation of the microscope images. The computational results presented have been achieved in part using the Vienna Scientific Cluster (VSC4).
Funding
Funding was provided by the Austrian Research Promotion Agency (FFG 877452 – INKplant).
Author information
Authors and Affiliations
Contributions
E.K. performed measurements and simulations and drafted the manuscript. M.K. and E.U. fabricated and prepared samples. L.K. and A.K. provided OCT scans. F.M. was supervising and reviewed the manuscript.
Corresponding author
Ethics declarations
Ethics approval and consent to participate
Not applicable.
Competing interests
The authors declare no competing interests.
Additional information
Publisher’s Note
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Appendices
Appendix A Multimaterial homogenization
In Eq. 1 is shown how the estimated compliance \(\textbf{C}^{est}\) of the multimaterial can be calculated from the compliance of the matrix material \(\textbf{C}_0\). The individual material phases r, s, have their respective volume fractions \(f_{r,s}\), compliances \(\textbf{C}_{r,s}\) and Hill tensors \(\textbf{P}_{0}^{r,s}\), where the geometrical information is inherited from the respective Eshelby tensor \(\textbf{S}_{0}^{r,s}\). The unit tensor is denoted by \(\mathbb {I}\).
If there is only one inclusion phase to consider, Eq. 1, can be simplified to Eq. 2 for the inclusion phase inc and the matrix phase 0. Mind, that the Hill tensor \(\textbf{P}\) considers the geometry of the inclusion, but the mechanical properties of the matrix material.
Estimation of the mechanical Eshelby tensor
Although the numerical simulations used the secondorder Eshelby tensor \(\textbf{S}\) in KelvinMandel notation, the mathematical derivation can be performed with the fourthorder Eshelby tensor \(\mathbb {S}_{MnAb}\) representation, as given in Eq. 3.
The secondorder tensor \(\varvec{K}_{MJ}\), dependent from the directions in the three dimensional space \(\vec {x}\), from Eq. 4 can be stated as in Eq. 7, assuming a compliance tensor \(\textbf{C}\) as in Eq. 6, where it is given in KelvinMandel notation. Together with the inclusions geometry, parameterized with the three dimensional space coordinates \(\vec {y}\), the surface integrals \(I_{abcd}\) can be solved, which are required to estimate the Eshelby tensor \(\mathbb {S}\).
The corresponding secondorder Hill tensor \(\textbf{P}\) can be calculated from the secondorder Eshelby tensor in KelvinMandel notation \(\textbf{S}\) and the material properties \(\textbf{C}\), according to Eq. 8.
Appendix B Sample materials
All used materials are commercially available, made for the use in the Connex3 Objet500 polyjet printer and products of Stratasys Ltd. (Minnesota, USA).

VeroClear is a transparent rigid material, also known as RGD810.

VeroPureWhite is an opaque rigid material, also known as RGD837.

TangoPlus is a transparent flexible material, also known as FLX930.

TangoBlackPlus is a flexible, opaque black material, also known as FLX980.

RGD8530DM is a “digital material” mixed in th printer during fabrication from VeroPureWhite and TangoBlackPlus.

RGD8430DM is a “digital material” mixed in the printer during fabrication from VeroPureWhite and TangoPlus.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Kornfellner, E., Königshofer, M., Krainz, L. et al. Measured and simulated mechanical properties of additively manufactured matrixinclusion multimaterials fabricated by material jetting. 3D Print Med 10, 4 (2024). https://doi.org/10.1186/s4120502300201y
Received:
Accepted:
Published:
DOI: https://doi.org/10.1186/s4120502300201y