Computational models for the simulation of the elastic and fracture properties of highly porous 3D‐printed hydroxyapatite scaffolds

Bone scaffolding is a promising approach for the treatment of critical‐size bone defects. Hydroxyapatite can be used to produce highly porous scaffolds as it mimics the mineralized part of bone tissue, but its intrinsic brittleness limits its usage. Among 3D printing techniques, vat photopolymerization allows for the best printing resolution for ceramic materials. In this study, we implemented a Computed micro‐Tomography based Finite Element Model of a hydroxyapatite porous scaffold fabricated by vat photopolymerization. We used the model in order to predict the elastic and fracture properties of the scaffold. From the stress–strain diagram of a simulated compression test, we computed the stiffness and the strength of the scaffolds. We found that three morphometric features substantially affect the crack pattern. In particular, the crack propagation is not only dependent on the trabecular thickness but also depends on the slenderness and orientation of the trabeculae with respect to the load. The results found in this study can be used for the design of ceramic scaffolds with heterogeneous pore distribution in order to tailor and predict the compressive strength.


| INTRODUCTION
Bone tissue remodeling is finely regulated by the human body through hormonal and mechanical stimuli. 1 Some pathologies and trauma may cause a critical size defect that cannot be restored directly by the body.An effective approach for the treatment of these conditions is the usage of scaffolds, which act as a template to promote tissue growth and regeneration. 2Although the autograft represents the gold standard for the scaffolding procedure, its shortage and complications due to the harvesting procedure make other solutions necessary. 3][6] Hydroxyapatite (HAP) is a ceramic material able to create a strong bond at the interface with the native bone since its composition is the same as the mineral component of the native bone tissue.Although synthetic HAP has been demonstrated to be thermodynamically stable in vivo, the presence of some ions, like in the native HAP, can enhance its degradation promoting the formation of new bone tissue. 7The main drawback that limits the usage of HAP is related to its intrinsic brittleness. 8dditive manufacturing (AM) techniques offer almost unbounded possibilities for the design of complex structures.0][11] Although the presence of pores reduces drastically the effective strength of the scaffold, high porosity is necessary in order to permit the growth and proliferation of cells, the fluid flow of nutrients, and waste removal.3][14] In a previous work by Fan et al, 15 the effect of the porosity on the strength and the Weibull modulus of HAP scaffolds have been deeply investigated.However, no other morphometric parameters have been considered for the fracture behavior.
The aim of this paper is to present computational models for assessing the strength and stiffness of VPP-derived HAP scaffolds designed for bone tissue engineering.Specifically, the paper will focus on the use of these models to evaluate the mechanical properties of the scaffold, including its compressive strength and elastic modulus.In particular, the nucleation and evolution of the cracks and their relationship with the morphometric features will be investigated in order to optimize the design of scaffolds with improved performance in bone tissue engineering applications.

| MATERIALS AND METHODS
The scaffolds analyzed in this work have been manufactured through VPP and thermally treated for overall 99 h, including multistage debinding and final sintering at 1300 C for 2 h 16 (Figure 1).
μCT scans using a Phoenix Nanotom S (Waygate Technologies/Baker Hughes Digital Solutions GmbH, Wunstorf, Germany) were then performed on the scaffolds: an acceleration voltage of 80 kV and a current of 120 μA were used for the acquisition, applying a 0:1 mm Cu X-ray filter.Rotation step size of 0:33 was used and four images were collected and integrated for each step.Acquisition time for each image was 1:5 s and the voxel size was 3:5 μm . 17In the present study, the μCT scans of two scaffolds printed on the basis of the same nominal geometry (i.e., same CAD model) have been used for morphological characterization and finite element simulations.The two scaffolds have a cylindrical shape with an approximate diameter of 5:5 mm and an approximate height of 8 mm.The inevitable changes in geometrical properties (such as trabecular thickness, pore size, and porosity) that will emerge between the two scaffolds have to be addressed in the manufacturing process since they were produced using identical geometrical models.The two samples will exhibit distinct mechanical responses as a result of these variations.
F I G U R E 1 Vat photopolymerization-fabricated samples: "green" sample obtained after printing and cleaning-(A) lateral view, (B) top view; (C) sintered hydroxyapatite scaffolds.
2.1 | From μCT images to hexahedral mesh μCT grayscale images with a pixel size of 3:5 μm have been subjected to a 2 pixels radius gaussian blur to reduce the noise and then resampled in the three cartesian directions to get a pixel size of 10:5 μm.Finally, the Otsu algorithm has been applied to get binary images. 18The 3D stack of binary images represents the input geometry of the simulation, where the voxel discretization corresponds to a hexahedral cartesian mesh.The final meshes are composed of approximately 43 and 37 millions of elements for the two scaffolds.

| Morphometric analyses
Because of the small distortion of the cylindrical shape of the samples caused by deformation during the sintering process, implying that a clear and well-defined cylindrical boundary cannot be identified, the porosity (ϕ) has been estimated using the following relationship where P w v represent the total number of white voxels (solid material), n z the number of slices of the 3D stack (767), that is, along the axial direction of the cylinders and b the total area of the circular cross-section.The value of b has been computed in three steps: (i) squeezing the 3D stack in a 2D image, adding up the value of the images constituting the 3D stack (Figure 2A); (ii) excluding the zero values (black pixels in Figure 2B) from the area of the rectangular matrix; and (iii) fitting the area occupied by the solid pixels with a circle with radius r.The area of the cross-section b is therefore the area enclosed by the fitted circle.
The trabecular thickness (Tb.Th), the trabecular spacing (Tb.Sp), and the ellipsoid factor (EF) of the scaffolds have been computed by using the BoneJ plugin of ImageJ. 19

| Macroscopic elastic properties
From both cylindrical samples, two cubes with edge of 300 pixels (3:15 mm) have been cropped from the cylindrical volume to determine the macroscopic elastic modulus of the porous material by applying homogenization approaches.The size of the cubes has been chosen in order to take the largest cube that fits into the cylindrical scaffolds.Therefore, only two cubes (one in the top region and one in the bottom region) fit into the entire cylindrical scaffold.Taking the largest cube fitting the entire sample will guarantee the representativeness of the cubic models to simulate the response of the material as a whole at the macroscopic level.Six units strain finite element simulations have been performed by using the solver ParOSol 20 to compute the effective homogenized elastic tensor components. 21The effective elastic moduli of the four cubic samples along the three cartesian directions have been compared for statistical analysis.Due to the small number of data, the assumption of non-normality has been done.Thereby, the U Mann-Whitney test has been used to find any significant statistical differences (significance value p = .05)between the two distributions of Young's moduli (mean elastic modulus in the xy plane (E xy mean =E 0 ) and E z =E 0 ) and between the two distributions of the shear moduli (mean shear modulus orthogonal to the xy plane (G xy z =E 0 ) and G xy =E 0 ).

| Macroscopic compressive strength
Uniaxial compression simulations have been performed on the whole cylindrical samples with the aim of simulating the post-elastic failure induced by the fracture propagation into the microstructures of the scaffolds.To this purpose, the ParOSol finite element solver was used within an iterative procedure able to simulate fracture propagation 22 ; the macroscopic strength was obtained as the maximum macroscopic compression stress before major failure of the samples occurred; this was identified as a significant drop in the macroscopic stress.On the bottom surface, the displacement component along the direction perpendicular to the plane was constrained; on the top surface the displacement component perpendicular to the plane with downward orientation (compression) was applied; the magnitude was progressively adjusted on the basis of the iterative procedure.On the lateral surface, two points were selected and displacement components belonging to the horizontal plane were constrained with the purpose of eliminating rigid body motions.The finite element model has been fed using the elastic constants E 0 ¼ 100 GPa 23 and assuming ν ¼ 0:3.The Drucker-Prager criterion has been used as strength criterion for local failure of the solid material, imposing the tensile-compressive mismatch by using σ T ¼ 100 MPa 23 and σ C ¼ 350 MPa 24 as intrinsic tensile strength and intrinsic compressive strength, respectively.The macroscopic stiffness of the whole cylindrical scaffold was also computed as the slope of the stress-strain diagram.

| RESULTS
Table 1 shows the major morphometric features of the two scaffolds.A trabecular thickness of approximately 300 μm and a trabecular spacing (i.e., pore size) of 800 μm were found; as reported by Abbasi et al 25 such spacing is adequate for bone cell proliferation. 26As a representative example, the trabecular spacing in the human femoral head is very similar to that in the scaffolds analyzed in this work, but variability in the human bone is more limited. 27The porosity of the two samples was 73% and 76%, respectively, which is also adequate for bone applications.
It is important to note that two of the key geometric design factors for a bone tissue engineering scaffold are the pore size and the trabecular thickness.The second stage of replacing bone tissue is where the pore size and interconnectivity become most important; they should ensure adequate fluid movement for nutrient transfer and waste elimination.This is anticipated to be accomplished by sufficiently large interconnected pores.As it must ensure proper stiffness and strength during the early stage loading, trabecular thickness is a factor playing a significant role in the first stage of the bone replacement process.In order to ensure adequate mechanical strength and simultaneously satisfy the technological requirements for its printability, an ideal trabecular thickness should be chosen.The porosity (i.e., the void-to-solid ratio) is the result of the simultaneous selection of pore size and trabecular thickness.
Figure 3A shows the spatial distribution of the trabecular thickness and, in panel (B) the corresponding histogram (blue) and the histogram of the trabecular spacing (orange).The x-axis represents trabecular thickness in micrometers (μm), and the y-axis represents the frequency of voxels belonging to a volume/section with a particular trabecular thickness or spacing.The blue histogram shows that the majority of the voxels in the image belong to volumes with The homogenized elastic moduli along the three cartesian directions and the shear stiffness in the three cartesian planes, as obtained through the homogenization procedure, are reported in Table 2; normalized values with respect to the elastic modulus of the bulk solid phase (E 0 ) are reported.For the sake of completeness, the porosity of each cubic domain is also reported.
The U Mann-Whitney tests show that there are no statistical differences between Young's moduli and shear moduli data in the three cartesian directions (p-value = .2for the former, p-value = .69for the latter).This indicates that  trabecular orientations are not statistically different among the three cartesian directions.On the basis of these statistical findings, the two scaffolds can be described using only three effective elastic constants (E, ν, and G), implying cubic symmetry.A quantitative measure of the anisotropy in a cubic material is given by the Zener ratio 28 (α z ): Averaging the four values of the Zener ratio related to cubic samples (α z ¼ 1:19 AE 0:02) reveals that the scaffolds are stiffer in shear than an isotropic equivalent structure (Figure 4).This result is consistent with the scaffold's architecture, which includes diagonal struts that act as braces.
The macroscopic stress versus macroscopic strain diagram of the two scaffolds under uniaxial compression as obtained by the numerical simulations is reported in Figure 5.The slope of the initial linear stage of the stress-strain response identifies the macroscopic elastic modulus along the loading direction which was 7:19 and 5:58 GPa, F I G U R E 4 Three-dimensional representation of the stiffness of the scaffold in the space.Since the scaffolds exhibit a cubic symmetry, only three material constants are needed (i.e., E=E 0 ¼ 0:063, G=E 0 ¼ 0:03, and ν ¼ 0:24).
F G U R E 5 Stress-strain diagram of the two scaffolds under uniaxial compression test.
respectively.The scaffold strength can be identified as the maximum stress achieved during the loading process which was σ C ¼ 3:34 MPa and σ C ¼ 2:53 MPa, respectively.
Figure 6 depicts the crack pattern after uniaxial compressive load.There are two types of fractures that have been observed: (i) cracks at the bifurcation and (ii) cracks at the interconnecting nodes.Fractures appear in both cases near the geometrical singularities of the scaffold, which act as stress amplifiers.
The scaffolds have pores of different sizes and shapes, which are not evenly distributed.As a result of this heterogeneous pore distribution, fractures occur in different parts of the scaffold rather than forming a damaged surface or narrow band.This finding is consistent with results by Wang et al 29 obtained for a Voronoi-based architecture.
Figure 7A reports two histograms.The first histogram (blue bars) shows the frequency distribution of voxels in the solid material separated into different classes based on the local trabecular thickness.This latter is obtained, at each solid voxel, as the maximum radius of the local sphere entirely contained into the solid material.
The second histogram (orange bars) shows the percentage of voxels belonging to the cracked volume for each class of trabecular thickness (i.e., the percentage of elements that underwent cracking and consequently removed from the computation).This histogram provides information about the distribution of cracks within the different trabecular thickness classes, which may be important for understanding the structural integrity and potential failure modes of the material.This latter histogram shows that fractures are less likely to occur in structures exhibiting a Tb.Th >500 μm.In particular, the distribution of the cracked volume is right-skewed, with a skewness of 10.9, whereas the trabecular thickness of the entire solid fraction has a skewness of 0.2, indicating a normal distribution.
The percentage of cracked volume in each thickness class is calculated by dividing the number of cracked elements in each thickness class by the solid volume in the same thickness class (Figure 7B).This histogram shows similar information to the left one but in terms of percentage.In each column, the percentage of elements (voxels) belonging to the particular thickness class that has been deleted (i.e., cracked) is reported.The highest number of deleted elements corresponds to thinner trabeculae, but there are other peaks at 210 and 315 μm.
Figure 8 shows the percentage of deleted elements in the portion of the solid with an EF greater than 0.66 (highly prolate ellipsoids, that is, thin rod-like structures).The angle formed by the load and the rod axis is used to group the bars.The graph demonstrates that trabecular structures with the axis aligned to the load in a range less than 20 are significantly less prone to fracture.This result is justified by the fact that trabecular structures misaligned with the loading direction are more subjected to bending and therefore to local tensile stress which eventually induces fracture and failure of the trabeculae.On the contrary, trabeculae aligned with the loading direction are more subjected to compressive stress and, therefore, less prone to fracture owing to the high compressive strength in comparison to the tensile strength.
F I G U R E 6 Rendering of the computational crack pattern under uniaxial compressive load.The crack pattern is highlighted in red.On the left-hand side, vertical cracks at the bifurcation; on the right-hand side, cracks perpendicular to the geometrical axis of the trabeculae.

| DISCUSSION
The aim of this paper is to use a computational model to determine the macroscopic compressive strength of highly porous scaffolds made of HAP created using the VPP technology.The highly porous nature of the scaffold makes it challenging to accurately determine its compressive strength experimentally.Therefore, a computational model can be useful in predicting the behavior of the scaffold under compressive loads.
The model takes into account the scaffold's geometry, pore size, and distribution, as well as the mechanical properties of the HAP material.By using the model to simulate the compressive loading of the scaffold, the macroscopic compressive strength can be determined.To this purpose, the geometry of two samples was determined through μCT.The μCT scans allowed us to determine the major morphometric parameters that affect the mechanical properties and a μCT-based finite element model was used to determine the macroscopic elastic and compressive strength properties.The finite element model was fed with the mechanical properties of the HAP as obtained from an experimental campaign suitably performed on miniaturized HAP samples created on purpose by the VPP technique. 23I G U R E 7 (A) Number of elements belonging to the solid fraction (in blue) and deleted elements (in orange) classified with respect to the Tb.Th; (B) ratio between the deleted elements and the total amount of solid fraction.
F I G U R E Bar plot of the percentage of deleted elements classified according to the angle between the geometrical axis of the rod and the load.
The computational models based on μCT scans have provided a macroscopic strength and macroscopic stiffness on the porous scaffold which is consistent with analytical predictive models available for porous materials. 30,31To this purpose, Figure 9 shows the normalized macroscopic elastic modulus E=E 0 in which E 0 is the intrinsic elastic modulus as in Reference 23 and the normalized compressive strength C T in which σ T is the intrinsic tensile strength as obtained in Reference 23.Comparison between the elastic modulus obtained from compression tests experimentally and numerically has to be done with caution.In particular, during experimental compression tests local fractures at the gripping area can decrease the value of the observed macroscopic stiffness.This makes the comparison between the experimental and numerical values of the elastic modulus not reliables.It has been already proven that nondestructive methods (i.e., acoustic methods) provide a good matching between the numeric and experimental macroscopic stiffness. 32,33igure 9 shows that although only a few data from numerical simulations (two scaffolds) are available, the computational results are in very good agreement with the analytical prediction.It is worth noticing that both scaffolds have been obtained from the same design, thereby the variations in morphometric and consequently mechanical features have to be addressed in the manufacturing process.The macroscopic strength of the two analyzed scaffolds was the result of a fracture pattern, shown in Figure 6, mostly generated by tensile stress; therefore, the macroscopic strength is weakly dependent on the assumed intrinsic compressive strength of the material.
The computational model slightly overestimates the compressive strength obtained by means of mechanical laboratory compression tests, that is, 1:60 AE 0:79 MPa on scaffolds with the same theoretical geometrical features. 16his is a key issue to consider when comparing computational and experimental results on porous HAP.While the computational model can provide a theoretical prediction of the compressive strength of a scaffold based on a specific geometry (as obtained by μCT data) and material properties, the experimental results may be affected by the actual manufacturing process and the resulting variations in the scaffold's geometry and microstructure.
In this case, the scaffolds used in the experimental tests were obtained through VPP technology with thermal treatment of the printed scaffolds, which can introduce geometrical deviations in trabecular thickness, shrinking, surface finishing or orientation of the microstructures, as well as variability in the intrinsic mechanical property of the solid material.These variations can affect the mechanical behavior of the scaffold under compressive loads and may result in different experimental results than those predicted by the computational model.Furthermore, imperfections of top and bottom surfaces of the real samples may heavily affect the experimental results.
Therefore, it is important to carefully consider the limitations of both the computational model and the experimental tests when interpreting and comparing the results.In particular, the variability introduced by the manufacturing process should be taken into account when evaluating the agreement between the computational and experimental results.
In order to have a better view of the comparison between experimental strength and the numerical prediction, the Weibull plot of the experiments reported in Figure 10 in Reference 16 is compared with pseudo-experimental results generated starting from the results of the computational model presented in this paper.In particular, 24 pseudoexperimental values of the compressive strength of the scaffolds were generated by extracting 12 random values for the F I G U R E 9 Comparison between the effective elastic modulus (square symbols) and the effective compressive strength (cross symbols) along the z-direction of the two samples and the polynomial fittings by Gibson and Ashby. 30,31ntrinsic strength σ T of HAP for each of the two scaffolds (there are a total of 24 random values); the random generation of intrinsic tensile strength was generated from a Weibull distribution with shape factor (m ¼ 2:6) and the scale factor (σ 0 ¼ 143) found by D'Andrea et al 23 on bulk 3D-printed miniaturized HAP samples.The computational prediction of the macroscopic compressive of each of the 12 intrinsic (for each of the two scaffolds) was obtained by scaling the macroscopic strength with the normalized strength shown in Figure 9.The 24 pseudo-experimental values are represented in Figure 10 along with the experimental values reported in Reference 16.The failure probability P F ¼ iÀ0: 5 24 for the pseudo-experimental values has been computed as defined by Munz and Fett.34 Figure 10 shows that the shape parameter for the pseudo-experimental strength is m ¼ 2:57 which is slightly higher than that experimentally found in Reference 16 (m ¼ 2:24).This mismatch is easily explained by the fact that the pseudo-generated results well catch variability in the intrinsic mechanical properties of the solid material, while variability in the architecture is not well represented as pseudo-experiments are obtained using only two μCT-based architectures.On the contrary, laboratory tests account for the variability of geometry and material simultaneously.Furthermore, the pseudo-experimental results exhibit a good overlap on the span of macroscopic strength between 1:2 and 3 MPa, but there is a rightward shift of the Weibull plot indicating a slight overestimation of the strength with respect to the experimental values obtained in samples with the same nominal porosity.The majority of the macroscopic strength values predicted by the FEM, as shown in this figure, overlap with the experimental values (between 1.5 and 2.7 MPa).This result confirms the validity of the proposed computational model despite the intrinsic limitation of this study due to the limited number of samples used in the numerical simulations.

| CONCLUSION
In this paper, the elastic and fracture characterization of VPP-derived HAP scaffolds through FEM have been reported.In particular, the role of the geometrical features of the structure in the cracking phenomenon has been investigated.Three main morphometric features are involved in the fracture process: the trabecular thickness, the ellipsoid factor and the orientation between the trabeculae and the load.In particular, fractures do not happen if the Tb.Th is larger than 500 μm and the probability of failure is lower in slender trabeculae aligned with the load.By choosing appropriate architectural features (such as porosity, trabecular thickness, and pore size), the proposed study's ultimate goal was to create a reliable numerical model that could be used in the design phase of the scaffold and be able to accurately predict the scaffold's properties within the variability caused by the manufacturing process.In order to increase the reliability and fidelity of the finished product, printing and manufacturing processes need to be improved.This goes beyond the scope of this particular study.The numerical model described in this work will be successfully integrated into an optimization design process for more dependable devices once this secondary, long-term goal is accomplished.
F I G U R E 1 0 Weibull distribution recomputed combining the experimental data from Baino et al 16 and the two computational strength data.

F
I G U R E 3 (A) spatial representation of the trabecular thickness, (B) histograms of the trabecular thickness (blue) and trabecular spacing (orange), (C) spatial representation of the ellipsoid factor and (D) histogram of the ellipsoid factor.T A B L E 2 Porosities homogenized elastic moduli.
Porosity, trabecular thickness, and the trabecular spacing of the two scaffolds.The trabecular thickness and the trabecular spacing are reported as mean value and standard deviation.and450μm, with a peak at approximately 350 μm.Figure3Cshows the spatial distribution of the ellipsoid factor and, in panel (D) the corresponding histogram.Plate-like trabeculae are represented by values close to À1, whereas rod-like trabeculae are represented by values close to 1.The histogram shows that the amount of solid material belonging to rod-like elements is greater than that belonging to plate-like structures, as expected in a foam-like scaffold.
T A B L E 1