- Email: [email protected]
Finite element modeling as a tool for predicting the fracture behavior of robocast scaffolds
Available online at www.sciencedirect.com
Acta Biomaterialia 4 (2008) 1715–1724 www.elsevier.com/locate/actabiomat
Finite element modeling as a tool for predicting the fracture behavior of robocast scaffolds Pedro Miranda *, Antonia Pajares, Fernando Guiberteau Departamento de Ingenierı´a Meca´nica, Energe´tica y de los Materiales, Universidad de Extremadura, Avda de Elvas s/n, 06071 Badajoz, Spain Received 10 January 2008; received in revised form 13 May 2008; accepted 22 May 2008 Available online 5 June 2008
Abstract The use of finite element modeling to calculate the stress fields in complex scaffold structures and thus predict their mechanical behavior during service (e.g., as load-bearing bone implants) is evaluated. The method is applied to identifying the fracture modes and estimating the strength of robocast hydroxyapatite and b-tricalcium phosphate scaffolds, consisting of a three-dimensional lattice of interpenetrating rods. The calculations are performed for three testing configurations: compression, tension and shear. Different testing orientations relative to the calcium phosphate rods are considered for each configuration. The predictions for the compressive configurations are compared to experimental data from uniaxial compression tests. Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Finite element analysis; Robocasting; Scaffolds; Strength; Fracture modes
1. Introduction Tissue engineering for bone regeneration requires the use of a porous scaffold to serve as a template for cell interactions and the formation of the extracellular matrix, as well as to provide structural support for the newly formed tissue [1]. These scaffolds should be able to withstand some degree of loading during their use in vivo, besides providing the required biological response. However, since porous scaffolds are inherently weak (especially those based on calcium phosphates) [2], it is necessary to optimize their mechanical response so that they can provide the appropriate support while bone regenerates and they are slowly resorbed [3]. Characterizing and predicting the mechanical behavior of these complex three-dimensional (3D) structures is therefore an essential task. Finite element modeling (FEM) has been shown to be capable of predicting the behavior of complex structures, such as multilayer systems [4–10], provided the mechanical
*
Corresponding author. Tel.: +34 924 28 9600; fax: +34 924 28 9601. E-mail address: [email protected] (P. Miranda).
properties of the materials comprising the structure are known. Consequently, the objective of the present work was to analyze the capabilities of finite element simulation as a technique for calculating the stress fields in complex scaffolds and thereby estimating their mechanical behavior under different loading configurations. In particular, FEM is applied to identify the fracture modes and estimate the strength of hydroxyapatite (HA) and b-tricalcium phosphate (b-TCP) scaffolds fabricated by robocasting. Robocasting, or direct-write assembly, is a solid freeform fabrication (SFF) technique that consists of the robotic deposition of highly concentrated colloidal suspensions (inks) capable of fully supporting their own weight during assembly [11–13]. Robocasting allows ceramic scaffolds to be built using water-based inks with minimal organic content (<1 wt.%) and without the need for a sacrificial support material or mold. Using this technique a network of semisolid ink rods is printed directly, hence the term direct-write assembly, by extrusion through a deposition nozzle to get the desired 3D structure. Fig. 1 shows a schematic diagram of the process. The prototype structure considered for this study consisted of a 3D tetragonal lattice of interpenetrating rods
1742-7061/$ - see front matter Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actbio.2008.05.020
1716
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 1. Scheme of the robocasting fabrication process. The ceramic scaffold is built layer by layer from a computer design. A three-axis robotic arm moves the injection syringe while pressing the ceramic ink through the conical deposition nozzle, immersed in an oil bath, to create a selfsupporting 3D network of ceramic rods.
– fabricated from either HA or b-TCP inks. The mechanical behavior of these scaffolds was simulated by FEM under three testing configurations: uniaxial compression, tension and shear. Different testing orientations relative to the rods comprising the structure were considered for each configuration. The most probable fracture mode was determined from the location of the maximum tensile stresses in the structure. This approach had been found to give excellent results in an earlier study of the compressive fracture modes of HA scaffolds [14]. The predicted compressive, tensile and shear strengths of the scaffolds for each testing orientation were estimated from the FEM results assuming a critical stress criterion. Finally, the predictions for the compressive configuration were compared to experimental data from uniaxial compression tests. 2. Experimental procedure Details of the fabrication of the HA and b-TCP inks and the corresponding robocast scaffolds have been presented elsewhere [2,15] and, as they are used here only for the experimental verification of the FEM results obtained in this work, only the relevant data will be included in this section. The 3D scaffolds, consisting of a lattice of ceramic rods, were constructed layer by layer via extrusion of the corre-
sponding inks through conical deposition nozzles, using a computer-controlled robotic deposition device (3D Inks, Stillwater, OK), as shown schematically in Fig. 1. The deposition was done in a non-wetting oil bath to prevent non-uniform drying during assembly. The samples were removed from the bath and dried in air at room temperature for 24 h, followed by a burn-out at 400 °C (ramp at 1 °C min 1 and a 1 h hold) to evaporate organics. They were finally sintered at 1300 °C (ramp at 3 °C min 1 and a 2 h hold). The resulting scaffolds were imaged (Fig. 2) using scanning electron microscopy (SEM) (S-4300SE/N, Hitachi, USA). The average rod diameter was measured to be d = 220 ± 10 lm and the average center-to-center spacing between rods in the printing plane s = 300 ± 10 lm, with a rod overlap between adjacent layers of about 50 ± 10 lm, yielding a vertical spacing between like layers of around h = 340 ± 10 lm (Fig. 3). The intrinsic mechanical properties of the individual calcium phosphate rods were evaluated to provide essential input parameters for the numerical simulations. Instrumented indentation (Nanotest, Micro Materials Ltd., Wrexham, UK) was used to determine the elastic moduli of the HA and b-TCP rods. Berkovich indentation tests were performed on polished sections (to 1 lm finish) perpendicular to the rod axis. Single indentations of about 5.5 lm depth and 40 lm side were placed at the center of the rods. This indent size is large enough compared to the grain sizes in both materials (3.2 ± 0.5 and 7.4 ± 0.7 lm for HA and b-TCP, respectively [2]) to provide meaningful information about the mechanical properties of the rods (and not of individual grains) but small enough to avoid a significant influence of the free surface of the rods – ASTM E 384 recommends a center-to-side spacing between Vickers indents of 2.5 times the diagonal. The inert fracture strength (i.e., the strength in the absence of chemical fatigue or slow crack growth) of the rods was estimated from high-load-resolution three-point bending tests (Electroforce 3200, Bose Corp., Eden Prairie, MN) performed on individual rods printed and sintered for the purpose. The tests were performed in air at a constant crosshead speed of 30 mm min–1. Since the fracture load for rods of 220 lm diameter was close to the sensitivity limit of the load cell, thicker rods of 360 lm diameter were used instead. Although the geometry and size of these isolated rods deviate from the real situation, these tests are the best possible means available to measure this material property (much more accurate than testing bulk samples of like materials, for example). Besides, in the opinion of the authors, the actual values of the rod strength should lay within the errors of these imperfect measurements. Finite element simulations were carried out using ABAQUS/StandardÒ software (Simulia, Providence, RI). The algorithm models an approximately cubic scaffold (2 2 2.09 mm) consisting of 12 alternating orthogonal layers of parallel HA or b-TCP rods of 220 lm diameter spaced 300 lm from center to center, with each layer overlapping the adjacent ones by 50 lm (Fig. 3). The FEM grid
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
1717
Fig. 2. SEM micrographs showing the morphology of HA (a, c) and b-TCP (b, d) scaffolds after sintering at 1300 °C for 2 h: (a, b) printing plane view and (c, d) cross-section view.
Fig. 3. Finite element grid used to simulate the scaffolds, consisting of 12 layers of rods of 220 lm diameter, d, spaced s = 300 lm from center-tocenter, and with a layer overlap of 50 lm that yields a like-layer spacing of h = 340 lm along the printing direction (direction 3). The dimensions of the elements at the rod surfaces are around 25 lm, but become larger towards the interior.
consists of more than 1.25 106 linear tetrahedral elements of about 25 lm in the vicinity of the external surfaces of the rods, but progressively increasing in size towards their interior. Isotropic elastic behavior is assumed for the whole
system, using the elastic moduli obtained from the indentation tests and a Poisson’s ratio of 0.28 [16] for both materials. The scaffold is placed between two parallel rigid plates (Fig. 4), already in contact with the structure from the beginning. The contact is considered frictionless for the compressive test simulations and infinitely rough and sticky for the tensile and shear test simulations. One of the rigid plates is fixed while the other moves under the action of a linearly increasing applied force (up to 500 N). For both the uniaxial tensile and compressive test simulations, two different orientations were considered, one with the load applied perpendicular to the printing plane (Fig. 4a and c) and the other with the load applied parallel to the rods (Fig. 4b and d), i.e., along the equivalent axis 1 or 2 in Fig. 3. For the shear test simulations the three different load orientations depicted in Fig. 4e–g were considered. The evolution of stresses with the applied stress – estimated by normalizing the applied load by the corresponding initial cross-sectional area of the simulated sample – was recorded during the simulation. The predicted strengths (compressive, tensile, or shear) of the scaffolds, i.e., the applied stress acting on the structure at failure, for each testing configuration were calculated assuming a critical stress criterion. First the greatest tensile stress acting at each point in the structure, rt, is determined from the FEM calculated stress tensor. Failure is then presumed to occur when the maximum of rt equals the inert fracture strength, rF, of the rods (as determined from the bending
1718
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 4. Scheme of the different testing configurations simulated in this study: compression along (a) the printing direction or (b) a rod axis; tension along (c) the printing direction or (d) a rod axis; shear (e) on the printing plane along a rod axis, or on a plane orthogonal to a rod axis along (f) the printing direction or (g) a rod axis. The rigid plates used in the FEM simulations to transfer the load to the structure are indicated. The bottom plate is held fixed while a load is applied to the top plate in the direction indicated by the arrows.
tests) – this implicitly assumes an opening fracture mode (mode I), which is that most commonly activated in brittle materials [17]. For the calculation of the maximum value of rt, the near-contact regions were ignored as they are not representative of the macroscopic behavior of the structure under each load configuration (i.e., uniaxial tension, compression or shear), but rather of the boundary conditions selected for the simulations. Furthermore, contact stresses will only produce localized damage [14] and would not be responsible for the ultimate failure of the structure except under high-loading-rate, e.g., impact, conditions. Experimental uniaxial compression tests were performed along the aforementioned two directions (Fig. 4a and b) on cubic blocks of about 2 mm side cut from the sintered HA and b-TCP scaffolds. The tests were carried out on a universal testing machine (AG-IS10kN, Shimadzu Corp., Kyoto, Japan) in air at a constant crosshead speed of 0.6 mm min 1 (the rate of 30 mm min 1 used in bending tests was impractical for the compression tests because it caused premature failure from contact damage and thus a slower rate was selected). The compressive strength of
the structure was estimated as the maximum applied load divided by the square external cross-section of the sample, and the results were compared with the FEM predictions. More than 10 samples were tested in each case in order to get statistically reliable values. 3. Results and discussion Table 1 summarizes the results of the characterization of the individual rods comprising the HA and b-TCP scaffolds. The elastic modulus, E, of each material was calculated from the effective modulus, E* = E/(1 m2), obtained from the indentation tests, assuming a value of Table 1 Intrinsic mechanical properties of HA and b-TCP rods
HA TCP a
E* (GPa)
ma
E (GPa)
rF (MPa)
83 ± 4 38 ± 8
0.28 0.28
82 ± 4 36 ± 7
68 ± 12 27 ± 9
From Ref. [16].
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
0.28 [16] for the Poisson’s ratio. As can be clearly appreciated, the HA rods exhibited more than twice the stiffness and strength of the b-TCP rods. These elastic properties were used as input for the FEM simulations of the different testing configurations depicted in Fig. 4. The FEM results confirmed that the stresses developed in the scaffolds for each configuration are obviously independent of the elastic modulus of the material considered, as shown for the case of uniaxial compression in Fig. 5, where the maximum of rt is represented vs. the applied stress for both HA and b-TCP. Consequently, the analysis that follows of the stress fields generated and the associated damage modes that are to be expected in each case is applicable to any scaffold with a similar geometry. Evidently, the strains do vary from one material to another, and the elastic modulus of the material comprising the individual rods will determine the effective modulus of the scaffold. Fig. 6 shows the contours of the greatest tensile stress, rt, developed in the scaffold during uniaxial compression at 250 N load (which is around the critical load for fracture in the experimental tests) for the two orientations considered (Fig. 4a and b). When compression is applied perpendicularly to the printing plane (Fig. 6a), i.e., along direction 3, rt presents a maximum at the centers of the unsupported segments – specifically at the top and bottom surfaces of those segments – of all the rods (although not so clearly visible in Fig. 6a for those aligned along direction 2). At those points, rt is directed along the corresponding rod axis (i.e., rt = r11 or rt = r22) and therefore, assuming an opening fracture mode (mode I), which, as mentioned above, is the most commonly activated in brittle materials [17], will
Fig. 5. Plot of maximum tensile stress in the structure (maximum of rt) vs. the applied stress (applied load normalized by initial cross-sectional area) for both HA and b-TCP scaffolds in the two compression directions considered: perpendicular to the printing plane (direction 3, Fig. 4a) and along the rod axes (direction 1 or 2, Fig. 4b). Note that the results are independent of the material.
1719
induce cracking perpendicular to the rod axes over the entire structure at those points. This type of cracking has indeed been observed in experimental compressive tests with this same orientation [14]. Of course, crack initiation does not necessarily occur at the exact location of the maximum tensile stresses because it also depends on the location of the largest flaws [18]. When the compressive load is applied along one of the rod axes (directions 1 or 2), the maximum tensile stresses are again located on the top and bottom surfaces of the horizontal rods (Fig. 6b), but now close to the joints with the vertical rods (see the inset). The orientation of these maximum tensile stresses is such that it will induce cracking (again, assuming mode I fracture) of the horizontal rods along those joints, effectively detaching them from the vertical rods that remain largely in compression. This prediction was also consistent with the curvature of the crack surfaces observed in the experimental tests [14]. The relative magnitude of the maximum stresses in the two orientations will be discussed in detail below. For the case of uniaxial tension tests, Fig. 7 shows the rt contours developed in the scaffold at 250 N load (for comparison with Fig. 6) for the two orientations considered (Fig. 4c and d). Note that, since stresses increase linearly with load (Fig. 5), contours at any other load can be calculated just by multiplying the stresses in the legend by the appropriate factor. When the loading is along direction 3 (Fig. 7a), and ignoring the values near the contacts, rt presents its maximum values at the joints between adjacent rod layers, as clearly shown in Fig. 7a for the rods aligned along direction 1 but also present in those aligned along direction 2. At a certain load, these maximum tensile stresses will induce fracture along those joints, separating each layer from its neighbors. When the tension is applied along one of the rod axes (directions 1 or 2), the tensile stresses are largely located on the vertical rods (Fig. 7b), with their maxima in the free segments of those rods close to (though not exactly at) the joints with the horizontal rods (see inset). The orientation of these maximum tensile stresses is largely along the rod axes, so that cracks originated by them will tend to break the vertical rods transversally. Fig. 8 shows analogous contour plots of rt at 250 N load for the three different shear configurations selected (Fig. 4e–g). For the sake of clarity, the direction of the applied load is included in the figures (recall that the bottom surface is kept fixed). Although the three configurations are significantly different, the stress contours generated exhibit many similarities. For all three orientations, the tensile stresses organize into bands of alternating intensity oriented at 45° with respect to the load axis. The maxima of these stresses, disregarding the near-contact field, are located at the joints between the different rod layers, the edges at the intersections acting as stress concentrators. Therefore, for the three configurations, cracking will occur at these locations and the different rod layers will tend to detach from each other. However, the actual cracking process will be slightly different in the case of the configuration
1720
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 6. FEM-generated stress contours corresponding to the highest tensile stress, rt, generated when a compressive load of 250 N is applied (a) orthogonal to the rods (direction 3) or (b) along a rod axis (direction 1, equivalent to direction 2).
in Fig. 8c, where adjacent layers tend to rotate relative to each other instead of simply one displacing the others, as in the cases of Fig. 8a and b. Concerning the relative intensities of the maximum stresses developed in each of the configurations analyzed, Fig. 9 shows the absolute maximum value of rt anywhere in the scaffold vs. the applied stress (i.e., applied load normalized by the corresponding scaffold’s external cross-section) for each testing configuration. This representation and the conclusions that derive from it are independent of the material considered and therefore have a universal validity for any scaffold with the same morphology. As can be clearly appreciated, there are very large differences in the maximum tensile stress levels that develop in the scaffolds depending on the type of load applied. As expected, compressive loads are the mildest, with stresses about three times lower than those that develop in the tension tests and more than an order of magnitude lower than those of the shear configurations. Of the two compression
configurations considered, the stresses are higher when the load is applied perpendicular to the printing plane (direction 3), but this situation is reversed under tensile loads, which yield higher stresses when applied along the rod axes (direction 1 or 2). The most deleterious configuration analyzed is that corresponding to shear load applied in planes perpendicular to the printing plane, along the direction of the rod axes (Fig. 4g). The torsion-like relative displacement of the rod layers generates a stress field around the joint that is significantly greater than in the two other shear configurations, which nonetheless both give stress levels more than twice as large as those of the tension configurations. An immediate, though anticipated, conclusion from these results is that brittle porous scaffolds fabricated by robocasting, regardless of the material they are made of, will exhibit far superior strength properties under compressive loads; tensile and shear stresses acting on them should be minimized, if not altogether suppressed, to improve their mechanical performance.
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
1721
Fig. 7. FEM-generated stress contours corresponding to the greatest tensile stress, rt, generated when a tensile load of 250 N is applied (a) orthogonal to the rods (direction 3) or (b) along a rod axis (direction 1, equivalent to direction 2).
From the plots in Fig. 9 it is simple to predict the mechanical strength of a scaffold for each testing configuration by using a critical stress criterion, provided the intrinsic strength of the constituent material is known. Indeed, these plots allow one to determine the applied stress at which the maximum tensile stress in the scaffold will equal the strength of the material, i.e., the structure strength under that type of load. Fig. 10 shows the strength values calculated following this procedure for the two materials considered in this study, HA and b-TCP, using the values of the inert strength, rF, from Table 1. The error bars represent standard deviations of the predictions due to uncertainties in the measurement of rF from the corresponding
three-point bending tests. As expected, the compressive strength values exhibited by both materials are significantly higher than their corresponding tensile and, especially, shear strengths. That the HA values are more than twice those corresponding to b-TCP is a simple reflection of their respective inert strengths. Finally, Fig. 11 shows the comparison between the predicted compressive strength and the experimental results obtained in the uniaxial compressive tests performed on the actual scaffolds for the two testing orientations considered (Fig. 4a and b). Again, the error bars represent the standard deviation of the data. The FEM predictions are seen to agree with the experimental compressive strengths
1722
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 8. FEM-generated stress contours corresponding to the greatest tensile stress, rt, generated when a shear load of 250 N is applied on the printing plane along a rod axis (direction 1, equivalent to direction 2) (a), or on a plane orthogonal to a rod axis (direction 1, equivalent to direction 2) along the printing direction (direction 3) (b) or along the perpendicular rod axis (direction 2, equivalent to direction 1) (c).
within uncertainties for both the HA and the b-TCP scaffolds. The apparently systematic slight (within the errors) overestimate of the compressive strength, especially notice-
able for b-TCP, could be attributed to the use of an inert strength value, while the actual experiments were performed at a low speed, where slow crack growth could be
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 9. Plot of maximum tensile stress in the scaffold (maximum of rt) vs. the applied stress (applied load normalized by initial cross-sectional area) for all the testing configurations analyzed: compression (two solid lines), tension (two dashed lines) and shear (three lines with dots). Legend notation is the same as used in Fig. 4.
1723
Fig. 11. Bar graph showing comparison between experimental measurements and FEM predictions of the compressive strength of HA and b-TCP scaffolds for the two compression orientations considered. Error bars represent standard deviations of the data.
fore, these results confirm the ability of FEM to predict the fracture behavior of robocast scaffolds both qualitatively (i.e., to predict their fracture modes [14]) and quantitatively (i.e., to calculate their strength), which was the main objective of the present work. 4. Conclusions and implications
Fig. 10. Bar graph showing HA and b-TCP scaffold strength values predicted by FEM for the different testing configurations, using the inert strength, rF, values of Table 1. Error bars represent standard deviations of predictions due to uncertainties in the measurement of rF. Label notation on the x-axis is the same as used in Fig. 4.
reducing the fracture strength of the scaffold. The TCP strength values used to calculate FEM estimates are less reliable than the corresponding HA values because, even after increasing the rod diameters, the critical loads for cracking lay close to the sensitivity threshold of the load cell used in the three-point bending test, which could explain the more significant deviation of FEM prediction from experimental values for this material. In sum, there-
Computer-aided fabrication methods, such as robocasting, are ideal for fabricating scaffolds with a designed pore structure that enables their biological and mechanical responses to be optimized to levels unattainable with other techniques. In particular, optimizing the mechanical performance of porous brittle scaffolds is critical since that is precisely what prevents their widespread application for regeneration of load-bearing bone defects. To perform such optimization, it is necessary to have a predictive tool as an aid in the design of the most favorable scaffold architecture. As the present results show, FEM could be that tool, as it enables the analysis and prediction of failure modes and the quantitative evaluation of important mechanical parameters, such as strength. In particular, the present study has successfully analyzed the stress fields and the expected fracture modes under seven different loading conditions (Fig. 4) of a robocast scaffold consisting of a tetragonal lattice of ceramic rods. From these results (Fig. 9), a key guideline for the design of brittle robocast scaffolds was derived: implant shape and location should be designed to minimize the tensile and, above all, shear stresses acting on the structure. This last type of load is extremely prejudicial for the mechanical performance of the scaffold as it generates intense tensile stresses that can lead to premature failure of the system.
1724
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fortunately, in most bone implant applications the loads will be largely compressive and therefore scaffold performance will not be so seriously jeopardized. For this particular type of load, another recommendation would be to preferentially orient the scaffolds so that the compression is parallel to one of the rod axes in order to minimize the tensile stresses generated in the structure. The FEM results were also used to predict the fracture strength of HA and b-TCP scaffolds fabricated by this technique under each type of load by applying a critical tensile stress criterion using their respective inert strengths, as measured in three-point bending tests, as the critical stress values. The compressive strengths predicted by the FEM simulations (for both materials in the two orientations considered) were successfully validated by comparison with experimental uniaxial compression test data, justifying the suitability of the present methodological approach for the optimization task at hand. Indeed, predictions obtained by FEM could allow one to determine how the different geometrical variables (rod thickness, spacing, interpenetration, relative angles, etc.) affect the mechanical behavior of the structure. Work in this direction is currently in progress. The results of those studies will make it possible to create an intelligent design of the scaffold’s initial computer-aided design model and thereby optimize its mechanical performance. Although we have focused here on the estimation of fracture modes and strength, nothing prevents this method from being applied to the optimization of other mechanical parameters, such as the elastic modulus of the structure. In sum, FEM is a powerful tool for predicting the mechanical performance of scaffolds fabricated by robocasting or other SFF techniques, and thus optimizing their mechanical performance through the intelligent design of their geometry, which is an unavoidable prerequisite for their use in loadbearing bone tissue engineering applications. Acknowledgements This work was supported by the Ministerio de Educacio´n y Ciencia (Spanish Government) and the Fondo Social Europeo (MAT2006-08720).
References [1] Karageorgiou V, Kaplan D. Porosity of 3D biomaterial scaffolds and osteogenesis. Biomaterials 2005;26:5474–91. [2] Miranda P, Pajares A, Saiz E, Tomsia AP, Guiberteau F. Mechanical properties of calcium phosphate scaffolds fabricated by robocasting. J Biomed Mater Res A 2008;85A:218–27. [3] Adachi T, Osako Y, Tanaka M, Hojo M, Hollister SJ. Framework for optimal design of porous scaffold microstructure by computational simulation of bone regeneration. Biomaterials 2006;27: 3964–72. [4] Zhao H, Miranda P, Lawn BR, Hu XZ. Cracking in ceramic/metal/ polymer trilayer systems. J Mater Res 2002;17:1102–11. [5] Chai H, Lawn BR. Cracking in brittle laminates from concentrated loads. Acta Mater 2002;50:2613–25. [6] Miranda P, Pajares A, Guiberteau F, Deng Y, Lawn BR. Designing damage-resistant brittle-coating structures: I. Bilayers. Acta Mater 2003;51:4347–56. [7] Miranda P, Pajares A, Guiberteau F, Deng Y, Zhao H, Lawn BR. Designing damage-resistant brittle-coating structures: II. Trilayers. Acta Mater 2003;51:4357–65. [8] Hsueh CH, Luttrell CR, Becher PF. Analyses of multilayered dental ceramics subjected to biaxial flexure tests. Dent Mater 2006;22:460–9. [9] Kim JH, Miranda P, Kim DK, Lawn BR. Effect of an adhesive interlayer on the fracture of a brittle coating on a supporting substrate. J Mater Res 2003;18:222–7. [10] Deng Y, Miranda P, Pajares A, Guiberteau F, Lawn BR. Fracture of ceramic/ceramic/polymer trilayers for biomechanical applications. J Biomed Mater Res A 2003;67A:828–33. [11] Cesarano III J, Segalman JR, Calvert P. Robocasting provides moldless fabrication from slurry deposition. Ceram Ind 1998;148:94–102. [12] Cesarano J, Calvert P. Freeforming objects with low-binder slurry. US Patent 6027326 (2000). [13] Smay JE, Cesarano J, Lewis JA. Colloidal inks for directed assembly of 3D periodic structures. Langmuir 2002;18:5429–37. [14] Miranda P, Pajares A, Saiz E, Tomsia AP, Guiberteau F. Fracture modes under uniaxial compression in hydroxyapatite scaffolds fabricated by robocasting. J Biomed Mater Res A 2007;83A:646–55. [15] Miranda P, Saiz E, Gryn K, Tomsia AP. Sintering and robocasting of beta-tricalcium phosphate scaffolds for orthopaedic applications. Acta Biomater 2006;2:457–66. [16] Grenoble DE, Dunn KL, Katz JL, Gilmore RS, Murty KL. Elastic properties of hard tissues and apatites. J Biomed Mater Res 1972;6:221–33. [17] Lawn BR. Fracture of brittle solids. Cambridge: Cambridge University Press; 1993. [18] Miranda P, Pajares A, Guiberteau F, Cumbrera FL, Lawn BR. Role of flaw statistics in contact fracture of brittle coatings. Acta Mater 2001;49:3719–26.
Acta Biomaterialia 4 (2008) 1715–1724 www.elsevier.com/locate/actabiomat
Finite element modeling as a tool for predicting the fracture behavior of robocast scaffolds Pedro Miranda *, Antonia Pajares, Fernando Guiberteau Departamento de Ingenierı´a Meca´nica, Energe´tica y de los Materiales, Universidad de Extremadura, Avda de Elvas s/n, 06071 Badajoz, Spain Received 10 January 2008; received in revised form 13 May 2008; accepted 22 May 2008 Available online 5 June 2008
Abstract The use of finite element modeling to calculate the stress fields in complex scaffold structures and thus predict their mechanical behavior during service (e.g., as load-bearing bone implants) is evaluated. The method is applied to identifying the fracture modes and estimating the strength of robocast hydroxyapatite and b-tricalcium phosphate scaffolds, consisting of a three-dimensional lattice of interpenetrating rods. The calculations are performed for three testing configurations: compression, tension and shear. Different testing orientations relative to the calcium phosphate rods are considered for each configuration. The predictions for the compressive configurations are compared to experimental data from uniaxial compression tests. Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. Keywords: Finite element analysis; Robocasting; Scaffolds; Strength; Fracture modes
1. Introduction Tissue engineering for bone regeneration requires the use of a porous scaffold to serve as a template for cell interactions and the formation of the extracellular matrix, as well as to provide structural support for the newly formed tissue [1]. These scaffolds should be able to withstand some degree of loading during their use in vivo, besides providing the required biological response. However, since porous scaffolds are inherently weak (especially those based on calcium phosphates) [2], it is necessary to optimize their mechanical response so that they can provide the appropriate support while bone regenerates and they are slowly resorbed [3]. Characterizing and predicting the mechanical behavior of these complex three-dimensional (3D) structures is therefore an essential task. Finite element modeling (FEM) has been shown to be capable of predicting the behavior of complex structures, such as multilayer systems [4–10], provided the mechanical
*
Corresponding author. Tel.: +34 924 28 9600; fax: +34 924 28 9601. E-mail address: [email protected] (P. Miranda).
properties of the materials comprising the structure are known. Consequently, the objective of the present work was to analyze the capabilities of finite element simulation as a technique for calculating the stress fields in complex scaffolds and thereby estimating their mechanical behavior under different loading configurations. In particular, FEM is applied to identify the fracture modes and estimate the strength of hydroxyapatite (HA) and b-tricalcium phosphate (b-TCP) scaffolds fabricated by robocasting. Robocasting, or direct-write assembly, is a solid freeform fabrication (SFF) technique that consists of the robotic deposition of highly concentrated colloidal suspensions (inks) capable of fully supporting their own weight during assembly [11–13]. Robocasting allows ceramic scaffolds to be built using water-based inks with minimal organic content (<1 wt.%) and without the need for a sacrificial support material or mold. Using this technique a network of semisolid ink rods is printed directly, hence the term direct-write assembly, by extrusion through a deposition nozzle to get the desired 3D structure. Fig. 1 shows a schematic diagram of the process. The prototype structure considered for this study consisted of a 3D tetragonal lattice of interpenetrating rods
1742-7061/$ - see front matter Ó 2008 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actbio.2008.05.020
1716
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 1. Scheme of the robocasting fabrication process. The ceramic scaffold is built layer by layer from a computer design. A three-axis robotic arm moves the injection syringe while pressing the ceramic ink through the conical deposition nozzle, immersed in an oil bath, to create a selfsupporting 3D network of ceramic rods.
– fabricated from either HA or b-TCP inks. The mechanical behavior of these scaffolds was simulated by FEM under three testing configurations: uniaxial compression, tension and shear. Different testing orientations relative to the rods comprising the structure were considered for each configuration. The most probable fracture mode was determined from the location of the maximum tensile stresses in the structure. This approach had been found to give excellent results in an earlier study of the compressive fracture modes of HA scaffolds [14]. The predicted compressive, tensile and shear strengths of the scaffolds for each testing orientation were estimated from the FEM results assuming a critical stress criterion. Finally, the predictions for the compressive configuration were compared to experimental data from uniaxial compression tests. 2. Experimental procedure Details of the fabrication of the HA and b-TCP inks and the corresponding robocast scaffolds have been presented elsewhere [2,15] and, as they are used here only for the experimental verification of the FEM results obtained in this work, only the relevant data will be included in this section. The 3D scaffolds, consisting of a lattice of ceramic rods, were constructed layer by layer via extrusion of the corre-
sponding inks through conical deposition nozzles, using a computer-controlled robotic deposition device (3D Inks, Stillwater, OK), as shown schematically in Fig. 1. The deposition was done in a non-wetting oil bath to prevent non-uniform drying during assembly. The samples were removed from the bath and dried in air at room temperature for 24 h, followed by a burn-out at 400 °C (ramp at 1 °C min 1 and a 1 h hold) to evaporate organics. They were finally sintered at 1300 °C (ramp at 3 °C min 1 and a 2 h hold). The resulting scaffolds were imaged (Fig. 2) using scanning electron microscopy (SEM) (S-4300SE/N, Hitachi, USA). The average rod diameter was measured to be d = 220 ± 10 lm and the average center-to-center spacing between rods in the printing plane s = 300 ± 10 lm, with a rod overlap between adjacent layers of about 50 ± 10 lm, yielding a vertical spacing between like layers of around h = 340 ± 10 lm (Fig. 3). The intrinsic mechanical properties of the individual calcium phosphate rods were evaluated to provide essential input parameters for the numerical simulations. Instrumented indentation (Nanotest, Micro Materials Ltd., Wrexham, UK) was used to determine the elastic moduli of the HA and b-TCP rods. Berkovich indentation tests were performed on polished sections (to 1 lm finish) perpendicular to the rod axis. Single indentations of about 5.5 lm depth and 40 lm side were placed at the center of the rods. This indent size is large enough compared to the grain sizes in both materials (3.2 ± 0.5 and 7.4 ± 0.7 lm for HA and b-TCP, respectively [2]) to provide meaningful information about the mechanical properties of the rods (and not of individual grains) but small enough to avoid a significant influence of the free surface of the rods – ASTM E 384 recommends a center-to-side spacing between Vickers indents of 2.5 times the diagonal. The inert fracture strength (i.e., the strength in the absence of chemical fatigue or slow crack growth) of the rods was estimated from high-load-resolution three-point bending tests (Electroforce 3200, Bose Corp., Eden Prairie, MN) performed on individual rods printed and sintered for the purpose. The tests were performed in air at a constant crosshead speed of 30 mm min–1. Since the fracture load for rods of 220 lm diameter was close to the sensitivity limit of the load cell, thicker rods of 360 lm diameter were used instead. Although the geometry and size of these isolated rods deviate from the real situation, these tests are the best possible means available to measure this material property (much more accurate than testing bulk samples of like materials, for example). Besides, in the opinion of the authors, the actual values of the rod strength should lay within the errors of these imperfect measurements. Finite element simulations were carried out using ABAQUS/StandardÒ software (Simulia, Providence, RI). The algorithm models an approximately cubic scaffold (2 2 2.09 mm) consisting of 12 alternating orthogonal layers of parallel HA or b-TCP rods of 220 lm diameter spaced 300 lm from center to center, with each layer overlapping the adjacent ones by 50 lm (Fig. 3). The FEM grid
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
1717
Fig. 2. SEM micrographs showing the morphology of HA (a, c) and b-TCP (b, d) scaffolds after sintering at 1300 °C for 2 h: (a, b) printing plane view and (c, d) cross-section view.
Fig. 3. Finite element grid used to simulate the scaffolds, consisting of 12 layers of rods of 220 lm diameter, d, spaced s = 300 lm from center-tocenter, and with a layer overlap of 50 lm that yields a like-layer spacing of h = 340 lm along the printing direction (direction 3). The dimensions of the elements at the rod surfaces are around 25 lm, but become larger towards the interior.
consists of more than 1.25 106 linear tetrahedral elements of about 25 lm in the vicinity of the external surfaces of the rods, but progressively increasing in size towards their interior. Isotropic elastic behavior is assumed for the whole
system, using the elastic moduli obtained from the indentation tests and a Poisson’s ratio of 0.28 [16] for both materials. The scaffold is placed between two parallel rigid plates (Fig. 4), already in contact with the structure from the beginning. The contact is considered frictionless for the compressive test simulations and infinitely rough and sticky for the tensile and shear test simulations. One of the rigid plates is fixed while the other moves under the action of a linearly increasing applied force (up to 500 N). For both the uniaxial tensile and compressive test simulations, two different orientations were considered, one with the load applied perpendicular to the printing plane (Fig. 4a and c) and the other with the load applied parallel to the rods (Fig. 4b and d), i.e., along the equivalent axis 1 or 2 in Fig. 3. For the shear test simulations the three different load orientations depicted in Fig. 4e–g were considered. The evolution of stresses with the applied stress – estimated by normalizing the applied load by the corresponding initial cross-sectional area of the simulated sample – was recorded during the simulation. The predicted strengths (compressive, tensile, or shear) of the scaffolds, i.e., the applied stress acting on the structure at failure, for each testing configuration were calculated assuming a critical stress criterion. First the greatest tensile stress acting at each point in the structure, rt, is determined from the FEM calculated stress tensor. Failure is then presumed to occur when the maximum of rt equals the inert fracture strength, rF, of the rods (as determined from the bending
1718
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 4. Scheme of the different testing configurations simulated in this study: compression along (a) the printing direction or (b) a rod axis; tension along (c) the printing direction or (d) a rod axis; shear (e) on the printing plane along a rod axis, or on a plane orthogonal to a rod axis along (f) the printing direction or (g) a rod axis. The rigid plates used in the FEM simulations to transfer the load to the structure are indicated. The bottom plate is held fixed while a load is applied to the top plate in the direction indicated by the arrows.
tests) – this implicitly assumes an opening fracture mode (mode I), which is that most commonly activated in brittle materials [17]. For the calculation of the maximum value of rt, the near-contact regions were ignored as they are not representative of the macroscopic behavior of the structure under each load configuration (i.e., uniaxial tension, compression or shear), but rather of the boundary conditions selected for the simulations. Furthermore, contact stresses will only produce localized damage [14] and would not be responsible for the ultimate failure of the structure except under high-loading-rate, e.g., impact, conditions. Experimental uniaxial compression tests were performed along the aforementioned two directions (Fig. 4a and b) on cubic blocks of about 2 mm side cut from the sintered HA and b-TCP scaffolds. The tests were carried out on a universal testing machine (AG-IS10kN, Shimadzu Corp., Kyoto, Japan) in air at a constant crosshead speed of 0.6 mm min 1 (the rate of 30 mm min 1 used in bending tests was impractical for the compression tests because it caused premature failure from contact damage and thus a slower rate was selected). The compressive strength of
the structure was estimated as the maximum applied load divided by the square external cross-section of the sample, and the results were compared with the FEM predictions. More than 10 samples were tested in each case in order to get statistically reliable values. 3. Results and discussion Table 1 summarizes the results of the characterization of the individual rods comprising the HA and b-TCP scaffolds. The elastic modulus, E, of each material was calculated from the effective modulus, E* = E/(1 m2), obtained from the indentation tests, assuming a value of Table 1 Intrinsic mechanical properties of HA and b-TCP rods
HA TCP a
E* (GPa)
ma
E (GPa)
rF (MPa)
83 ± 4 38 ± 8
0.28 0.28
82 ± 4 36 ± 7
68 ± 12 27 ± 9
From Ref. [16].
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
0.28 [16] for the Poisson’s ratio. As can be clearly appreciated, the HA rods exhibited more than twice the stiffness and strength of the b-TCP rods. These elastic properties were used as input for the FEM simulations of the different testing configurations depicted in Fig. 4. The FEM results confirmed that the stresses developed in the scaffolds for each configuration are obviously independent of the elastic modulus of the material considered, as shown for the case of uniaxial compression in Fig. 5, where the maximum of rt is represented vs. the applied stress for both HA and b-TCP. Consequently, the analysis that follows of the stress fields generated and the associated damage modes that are to be expected in each case is applicable to any scaffold with a similar geometry. Evidently, the strains do vary from one material to another, and the elastic modulus of the material comprising the individual rods will determine the effective modulus of the scaffold. Fig. 6 shows the contours of the greatest tensile stress, rt, developed in the scaffold during uniaxial compression at 250 N load (which is around the critical load for fracture in the experimental tests) for the two orientations considered (Fig. 4a and b). When compression is applied perpendicularly to the printing plane (Fig. 6a), i.e., along direction 3, rt presents a maximum at the centers of the unsupported segments – specifically at the top and bottom surfaces of those segments – of all the rods (although not so clearly visible in Fig. 6a for those aligned along direction 2). At those points, rt is directed along the corresponding rod axis (i.e., rt = r11 or rt = r22) and therefore, assuming an opening fracture mode (mode I), which, as mentioned above, is the most commonly activated in brittle materials [17], will
Fig. 5. Plot of maximum tensile stress in the structure (maximum of rt) vs. the applied stress (applied load normalized by initial cross-sectional area) for both HA and b-TCP scaffolds in the two compression directions considered: perpendicular to the printing plane (direction 3, Fig. 4a) and along the rod axes (direction 1 or 2, Fig. 4b). Note that the results are independent of the material.
1719
induce cracking perpendicular to the rod axes over the entire structure at those points. This type of cracking has indeed been observed in experimental compressive tests with this same orientation [14]. Of course, crack initiation does not necessarily occur at the exact location of the maximum tensile stresses because it also depends on the location of the largest flaws [18]. When the compressive load is applied along one of the rod axes (directions 1 or 2), the maximum tensile stresses are again located on the top and bottom surfaces of the horizontal rods (Fig. 6b), but now close to the joints with the vertical rods (see the inset). The orientation of these maximum tensile stresses is such that it will induce cracking (again, assuming mode I fracture) of the horizontal rods along those joints, effectively detaching them from the vertical rods that remain largely in compression. This prediction was also consistent with the curvature of the crack surfaces observed in the experimental tests [14]. The relative magnitude of the maximum stresses in the two orientations will be discussed in detail below. For the case of uniaxial tension tests, Fig. 7 shows the rt contours developed in the scaffold at 250 N load (for comparison with Fig. 6) for the two orientations considered (Fig. 4c and d). Note that, since stresses increase linearly with load (Fig. 5), contours at any other load can be calculated just by multiplying the stresses in the legend by the appropriate factor. When the loading is along direction 3 (Fig. 7a), and ignoring the values near the contacts, rt presents its maximum values at the joints between adjacent rod layers, as clearly shown in Fig. 7a for the rods aligned along direction 1 but also present in those aligned along direction 2. At a certain load, these maximum tensile stresses will induce fracture along those joints, separating each layer from its neighbors. When the tension is applied along one of the rod axes (directions 1 or 2), the tensile stresses are largely located on the vertical rods (Fig. 7b), with their maxima in the free segments of those rods close to (though not exactly at) the joints with the horizontal rods (see inset). The orientation of these maximum tensile stresses is largely along the rod axes, so that cracks originated by them will tend to break the vertical rods transversally. Fig. 8 shows analogous contour plots of rt at 250 N load for the three different shear configurations selected (Fig. 4e–g). For the sake of clarity, the direction of the applied load is included in the figures (recall that the bottom surface is kept fixed). Although the three configurations are significantly different, the stress contours generated exhibit many similarities. For all three orientations, the tensile stresses organize into bands of alternating intensity oriented at 45° with respect to the load axis. The maxima of these stresses, disregarding the near-contact field, are located at the joints between the different rod layers, the edges at the intersections acting as stress concentrators. Therefore, for the three configurations, cracking will occur at these locations and the different rod layers will tend to detach from each other. However, the actual cracking process will be slightly different in the case of the configuration
1720
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 6. FEM-generated stress contours corresponding to the highest tensile stress, rt, generated when a compressive load of 250 N is applied (a) orthogonal to the rods (direction 3) or (b) along a rod axis (direction 1, equivalent to direction 2).
in Fig. 8c, where adjacent layers tend to rotate relative to each other instead of simply one displacing the others, as in the cases of Fig. 8a and b. Concerning the relative intensities of the maximum stresses developed in each of the configurations analyzed, Fig. 9 shows the absolute maximum value of rt anywhere in the scaffold vs. the applied stress (i.e., applied load normalized by the corresponding scaffold’s external cross-section) for each testing configuration. This representation and the conclusions that derive from it are independent of the material considered and therefore have a universal validity for any scaffold with the same morphology. As can be clearly appreciated, there are very large differences in the maximum tensile stress levels that develop in the scaffolds depending on the type of load applied. As expected, compressive loads are the mildest, with stresses about three times lower than those that develop in the tension tests and more than an order of magnitude lower than those of the shear configurations. Of the two compression
configurations considered, the stresses are higher when the load is applied perpendicular to the printing plane (direction 3), but this situation is reversed under tensile loads, which yield higher stresses when applied along the rod axes (direction 1 or 2). The most deleterious configuration analyzed is that corresponding to shear load applied in planes perpendicular to the printing plane, along the direction of the rod axes (Fig. 4g). The torsion-like relative displacement of the rod layers generates a stress field around the joint that is significantly greater than in the two other shear configurations, which nonetheless both give stress levels more than twice as large as those of the tension configurations. An immediate, though anticipated, conclusion from these results is that brittle porous scaffolds fabricated by robocasting, regardless of the material they are made of, will exhibit far superior strength properties under compressive loads; tensile and shear stresses acting on them should be minimized, if not altogether suppressed, to improve their mechanical performance.
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
1721
Fig. 7. FEM-generated stress contours corresponding to the greatest tensile stress, rt, generated when a tensile load of 250 N is applied (a) orthogonal to the rods (direction 3) or (b) along a rod axis (direction 1, equivalent to direction 2).
From the plots in Fig. 9 it is simple to predict the mechanical strength of a scaffold for each testing configuration by using a critical stress criterion, provided the intrinsic strength of the constituent material is known. Indeed, these plots allow one to determine the applied stress at which the maximum tensile stress in the scaffold will equal the strength of the material, i.e., the structure strength under that type of load. Fig. 10 shows the strength values calculated following this procedure for the two materials considered in this study, HA and b-TCP, using the values of the inert strength, rF, from Table 1. The error bars represent standard deviations of the predictions due to uncertainties in the measurement of rF from the corresponding
three-point bending tests. As expected, the compressive strength values exhibited by both materials are significantly higher than their corresponding tensile and, especially, shear strengths. That the HA values are more than twice those corresponding to b-TCP is a simple reflection of their respective inert strengths. Finally, Fig. 11 shows the comparison between the predicted compressive strength and the experimental results obtained in the uniaxial compressive tests performed on the actual scaffolds for the two testing orientations considered (Fig. 4a and b). Again, the error bars represent the standard deviation of the data. The FEM predictions are seen to agree with the experimental compressive strengths
1722
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 8. FEM-generated stress contours corresponding to the greatest tensile stress, rt, generated when a shear load of 250 N is applied on the printing plane along a rod axis (direction 1, equivalent to direction 2) (a), or on a plane orthogonal to a rod axis (direction 1, equivalent to direction 2) along the printing direction (direction 3) (b) or along the perpendicular rod axis (direction 2, equivalent to direction 1) (c).
within uncertainties for both the HA and the b-TCP scaffolds. The apparently systematic slight (within the errors) overestimate of the compressive strength, especially notice-
able for b-TCP, could be attributed to the use of an inert strength value, while the actual experiments were performed at a low speed, where slow crack growth could be
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fig. 9. Plot of maximum tensile stress in the scaffold (maximum of rt) vs. the applied stress (applied load normalized by initial cross-sectional area) for all the testing configurations analyzed: compression (two solid lines), tension (two dashed lines) and shear (three lines with dots). Legend notation is the same as used in Fig. 4.
1723
Fig. 11. Bar graph showing comparison between experimental measurements and FEM predictions of the compressive strength of HA and b-TCP scaffolds for the two compression orientations considered. Error bars represent standard deviations of the data.
fore, these results confirm the ability of FEM to predict the fracture behavior of robocast scaffolds both qualitatively (i.e., to predict their fracture modes [14]) and quantitatively (i.e., to calculate their strength), which was the main objective of the present work. 4. Conclusions and implications
Fig. 10. Bar graph showing HA and b-TCP scaffold strength values predicted by FEM for the different testing configurations, using the inert strength, rF, values of Table 1. Error bars represent standard deviations of predictions due to uncertainties in the measurement of rF. Label notation on the x-axis is the same as used in Fig. 4.
reducing the fracture strength of the scaffold. The TCP strength values used to calculate FEM estimates are less reliable than the corresponding HA values because, even after increasing the rod diameters, the critical loads for cracking lay close to the sensitivity threshold of the load cell used in the three-point bending test, which could explain the more significant deviation of FEM prediction from experimental values for this material. In sum, there-
Computer-aided fabrication methods, such as robocasting, are ideal for fabricating scaffolds with a designed pore structure that enables their biological and mechanical responses to be optimized to levels unattainable with other techniques. In particular, optimizing the mechanical performance of porous brittle scaffolds is critical since that is precisely what prevents their widespread application for regeneration of load-bearing bone defects. To perform such optimization, it is necessary to have a predictive tool as an aid in the design of the most favorable scaffold architecture. As the present results show, FEM could be that tool, as it enables the analysis and prediction of failure modes and the quantitative evaluation of important mechanical parameters, such as strength. In particular, the present study has successfully analyzed the stress fields and the expected fracture modes under seven different loading conditions (Fig. 4) of a robocast scaffold consisting of a tetragonal lattice of ceramic rods. From these results (Fig. 9), a key guideline for the design of brittle robocast scaffolds was derived: implant shape and location should be designed to minimize the tensile and, above all, shear stresses acting on the structure. This last type of load is extremely prejudicial for the mechanical performance of the scaffold as it generates intense tensile stresses that can lead to premature failure of the system.
1724
P. Miranda et al. / Acta Biomaterialia 4 (2008) 1715–1724
Fortunately, in most bone implant applications the loads will be largely compressive and therefore scaffold performance will not be so seriously jeopardized. For this particular type of load, another recommendation would be to preferentially orient the scaffolds so that the compression is parallel to one of the rod axes in order to minimize the tensile stresses generated in the structure. The FEM results were also used to predict the fracture strength of HA and b-TCP scaffolds fabricated by this technique under each type of load by applying a critical tensile stress criterion using their respective inert strengths, as measured in three-point bending tests, as the critical stress values. The compressive strengths predicted by the FEM simulations (for both materials in the two orientations considered) were successfully validated by comparison with experimental uniaxial compression test data, justifying the suitability of the present methodological approach for the optimization task at hand. Indeed, predictions obtained by FEM could allow one to determine how the different geometrical variables (rod thickness, spacing, interpenetration, relative angles, etc.) affect the mechanical behavior of the structure. Work in this direction is currently in progress. The results of those studies will make it possible to create an intelligent design of the scaffold’s initial computer-aided design model and thereby optimize its mechanical performance. Although we have focused here on the estimation of fracture modes and strength, nothing prevents this method from being applied to the optimization of other mechanical parameters, such as the elastic modulus of the structure. In sum, FEM is a powerful tool for predicting the mechanical performance of scaffolds fabricated by robocasting or other SFF techniques, and thus optimizing their mechanical performance through the intelligent design of their geometry, which is an unavoidable prerequisite for their use in loadbearing bone tissue engineering applications. Acknowledgements This work was supported by the Ministerio de Educacio´n y Ciencia (Spanish Government) and the Fondo Social Europeo (MAT2006-08720).
References [1] Karageorgiou V, Kaplan D. Porosity of 3D biomaterial scaffolds and osteogenesis. Biomaterials 2005;26:5474–91. [2] Miranda P, Pajares A, Saiz E, Tomsia AP, Guiberteau F. Mechanical properties of calcium phosphate scaffolds fabricated by robocasting. J Biomed Mater Res A 2008;85A:218–27. [3] Adachi T, Osako Y, Tanaka M, Hojo M, Hollister SJ. Framework for optimal design of porous scaffold microstructure by computational simulation of bone regeneration. Biomaterials 2006;27: 3964–72. [4] Zhao H, Miranda P, Lawn BR, Hu XZ. Cracking in ceramic/metal/ polymer trilayer systems. J Mater Res 2002;17:1102–11. [5] Chai H, Lawn BR. Cracking in brittle laminates from concentrated loads. Acta Mater 2002;50:2613–25. [6] Miranda P, Pajares A, Guiberteau F, Deng Y, Lawn BR. Designing damage-resistant brittle-coating structures: I. Bilayers. Acta Mater 2003;51:4347–56. [7] Miranda P, Pajares A, Guiberteau F, Deng Y, Zhao H, Lawn BR. Designing damage-resistant brittle-coating structures: II. Trilayers. Acta Mater 2003;51:4357–65. [8] Hsueh CH, Luttrell CR, Becher PF. Analyses of multilayered dental ceramics subjected to biaxial flexure tests. Dent Mater 2006;22:460–9. [9] Kim JH, Miranda P, Kim DK, Lawn BR. Effect of an adhesive interlayer on the fracture of a brittle coating on a supporting substrate. J Mater Res 2003;18:222–7. [10] Deng Y, Miranda P, Pajares A, Guiberteau F, Lawn BR. Fracture of ceramic/ceramic/polymer trilayers for biomechanical applications. J Biomed Mater Res A 2003;67A:828–33. [11] Cesarano III J, Segalman JR, Calvert P. Robocasting provides moldless fabrication from slurry deposition. Ceram Ind 1998;148:94–102. [12] Cesarano J, Calvert P. Freeforming objects with low-binder slurry. US Patent 6027326 (2000). [13] Smay JE, Cesarano J, Lewis JA. Colloidal inks for directed assembly of 3D periodic structures. Langmuir 2002;18:5429–37. [14] Miranda P, Pajares A, Saiz E, Tomsia AP, Guiberteau F. Fracture modes under uniaxial compression in hydroxyapatite scaffolds fabricated by robocasting. J Biomed Mater Res A 2007;83A:646–55. [15] Miranda P, Saiz E, Gryn K, Tomsia AP. Sintering and robocasting of beta-tricalcium phosphate scaffolds for orthopaedic applications. Acta Biomater 2006;2:457–66. [16] Grenoble DE, Dunn KL, Katz JL, Gilmore RS, Murty KL. Elastic properties of hard tissues and apatites. J Biomed Mater Res 1972;6:221–33. [17] Lawn BR. Fracture of brittle solids. Cambridge: Cambridge University Press; 1993. [18] Miranda P, Pajares A, Guiberteau F, Cumbrera FL, Lawn BR. Role of flaw statistics in contact fracture of brittle coatings. Acta Mater 2001;49:3719–26.