Next Article in Journal
An Extension of Explicit Coupling for Fluid–Structure Interaction Problems
Next Article in Special Issue
Obtaining Expressions for Time-Dependent Functions That Describe the Unsteady Properties of Swirling Jets of Viscous Fluid
Previous Article in Journal
On State Occupancies, First Passage Times and Duration in Non-Homogeneous Semi-Markov Chains
Previous Article in Special Issue
Numerical Assessment of the Structural Effects of Relative Sliding between Tissues in a Finite Element Model of the Foot
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Size Effects in Finite Element Modelling of 3D Printed Bone Scaffolds Using Hydroxyapatite PEOT/PBT Composites

1
TECNALIA, Basque Research and Technology Alliance (BRTA), Mikeletegi Pasealekua 2, 20009 Donostia-San Sebastián, Spain
2
Department of Engineering, Public University of Navarra, 31006 Pamplona, Spain
3
Complex Tissue Regeneration Department, MERLN Institute for Technology-Inspired Regenerative Medicine, University of Maastricht, 6211 LK Maastricht, The Netherlands
4
Department of Physics and Astronomy, Padova University, Via Marzolo 8, 35131 Padova, Italy
5
Materials Physics Center (MPC), Centro de Fisica de Materiales (CSIC, UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain
*
Authors to whom correspondence should be addressed.
Submission received: 2 June 2021 / Revised: 6 July 2021 / Accepted: 7 July 2021 / Published: 24 July 2021
(This article belongs to the Special Issue Numerical Simulation in Biomechanics and Biomedical Engineering)

Abstract

:
Additive manufacturing (AM) of scaffolds enables the fabrication of customized patient-specific implants for tissue regeneration. Scaffold customization does not involve only the macroscale shape of the final implant, but also their microscopic pore geometry and material properties, which are dependent on optimizable topology. A good match between the experimental data of AM scaffolds and the models is obtained when there is just a few millimetres at least in one direction. Here, we describe a methodology to perform finite element modelling on AM scaffolds for bone tissue regeneration with clinically relevant dimensions (i.e., volume > 1 cm3). The simulation used an equivalent cubic eight node finite elements mesh, and the materials properties were derived both empirically and numerically, from bulk material direct testing and simulated tests on scaffolds. The experimental validation was performed using poly(ethylene oxide terephthalate)-poly(butylene terephthalate) (PEOT/PBT) copolymers and 45 wt% nano hydroxyapatite fillers composites. By applying this methodology on three separate scaffold architectures with volumes larger than 1 cm3, the simulations overestimated the scaffold performance, resulting in 150–290% stiffer than average values obtained in the validation tests. The results mismatch highlighted the relevance of the lack of printing accuracy that is characteristic of the additive manufacturing process. Accordingly, a sensitivity analysis was performed on nine detected uncertainty sources, studying their influence. After the definition of acceptable execution tolerances and reliability levels, a design factor was defined to calibrate the methodology under expectable and conservative scenarios.

1. Introduction

In order to manufacture patient-specific bone implants for tissue regeneration (Figure 1), it is important to keep in mind that any bone is a composite structure subjected to constant evolution normally composed of a high percentage of inorganic composition and a smaller organic fraction. The inorganic part is present intrafibrillarly in the bone and occupies up to 40% by volume filling with various degrees of space around the mineralized fibres, forming pore networks. Hence, human bones can be composed of high amounts of inorganics by volume, typically ranging from 50–60% up to 80–90% for some highly mineralized tissues, but always keeping some space for organics [1].
Consequently, with these boundary conditions, the engineering of patient-specific bone implants for tissue regeneration follows several approaches [2,3]. One of the most promising consists of the additive manufacturing of a temporary structure also known as a “scaffold”, suitable for cell growth, with required porosity to stimulate osteogenesis, as stated by Jakus et al. [4], Karageorgiou and Kaplan [5] and Hutmacher [6]. These temporary scaffolds must meet several requirements, such as biocompatibility, biodegradability and suitable mechanical properties [2,3,7,8,9,10,11,12,13], but also printability [14,15]. The mechanical properties are of great importance, since the scaffolds should ideally have properties matching those of the surrounding tissue, ensuring its suitability to bear mechanical loading similar to the original tissue. The mechanical properties are mainly dictated by the bulk material properties. However, the geometry and porosity of the fabricated scaffolds also play an important role in the final mechanical properties of the implant [16,17].
There are diverse materials being investigated, spanning from synthetic and natural polymers, ceramics and combinations of these, with the research focused on finding an optimal formulation for a successful stimulation of bone healing [2,3], enabling the improvement of implants suitable for each patient condition in the future. Thereby, additive manufacturing of a bone scaffold requires a compatible biocomposite with reinforcing inorganic fillers to mimic natural bone composition and structure, while improving its mechanical behaviour to resist service life loadings [18,19]. Accordingly, there are many studies regarding material alternatives [20], such as silicate nanocomposites [21,22], graphene nanocomposites [23] and the extensive work on hydroxyapatite or nano-apatite [24,25,26,27,28]. Additionally, other studies focused on the mechanical optimization and performance of such material alternatives in the presence of micro and nano hydroxyapatite [29], or on bone regeneration and thermomechanical properties [30,31].
In this manuscript, a generalized procedure is developed for the numerical modelling of additive manufactured scaffolds for bone regeneration with clinically relevant dimensions (i.e., volume > 1 cm3). The procedure is applicable to scaffolds made of any previously mechanically characterized material, and this procedure is then specifically validated with poly(ethylene oxide terephthalate)/poly(butylene terephthalate) (PEOT/PBT) with nano-hydroxyapatite fillers. This generalized procedure allows the scaffolds to be designed by function, fixing at first the requirements of geometry and the stresses to be supported and then selecting the optimal filler amount and porosity. A stiffer scaffold, for example, would require a greater amount of filler and less pores for stress to reduce the concentration spots. Consequently, this leads to design gradients and patterns for a single bone scaffold, which can be achieved using established topology optimization techniques.
Figure 1. 3D printed bone scaffolds for tissue regeneration: (a) long bone defect-shaped scaffold printed with PEOT/PBT alone; (b) Scaffolds presenting longitudinal gradient patterns of plasma polymerized film deposition, visualized using methylene blue staining [32].
Figure 1. 3D printed bone scaffolds for tissue regeneration: (a) long bone defect-shaped scaffold printed with PEOT/PBT alone; (b) Scaffolds presenting longitudinal gradient patterns of plasma polymerized film deposition, visualized using methylene blue staining [32].
Mathematics 09 01746 g001
The mechanical properties of the bone scaffolds are a key aspect for the development of PEOT/PBT composites and a fast-complete characterization of these bulk materials can be carried out using monotonic tensile and compression tests according to standards. Nevertheless, the equivalent testing of non-standardized bone scaffolds made of such materials is not as easy, for the scaffolds need to be previously manufactured with different geometries and architectures, which slow down their rapid characterization. Therefore, a modelling procedure suitable to predict the mechanical performance of bone scaffolds with a certain pore architecture and bulk material mechanical properties has been developed.
Finally, a set of scaffolds has been produced and mechanically characterised as a way to confirm the suitability of the developed modelling protocol and its validation. Nevertheless, the 3D printing of clinically relevant scaffolds with the required micro-scale precision is a challenging task and the manufacturing process generally produces unintended imperfections causing uncertainty because of result scattering. Hence, a sensitivity analysis and statistical study have been conducted in order to derive some design factors taking into account these issues.
The final output is a finite element modelling (FEM) procedure for 3D-printed bone scaffolds able to forecast their mechanical behaviour depending on the bulk material properties. This will enable the optimization of the material development process and perform topology optimization of such scaffolds for bone tissues. This is not a topology optimization necessarily based on evolutionary structural optimization (ESO) algorithms, removing material where it is not working at certain stress level; nor it is based on additive evolutionary structural optimization (AESO) algorithms, simply adding materials where the stresses are above certain stress level, since the global shape of the bone implant needs to be kept “custom made”, performing an external geometry for a specific patient. It is a topology optimization based on the scaffold properties’ distribution within such a patient custom made geometry, i.e., manufacturing implants with stronger bone scaffolds at stress concentration “hot spots” by varying layer height, strand distance, fibre diameter (changing extrusion nozzle or printing speed), bulk material properties, plasma treatments, etc.
There are several studies regarding the FEM of 3D bone scaffolds for tissue regeneration. For instance, there are several bone scaffold aspects being currently analysed in terms of scaffold design [33], mechanical behaviour during failure [34], or covering the bone tissue regeneration and the progressive change in properties [35,36]. Nevertheless, the scope of such studies is focused on the mechanical behaviour of the 3D bone scaffold itself. Yet, there is little research on the holistic methodology connecting the bulk material properties with the scaffold mechanical behaviour within a 3D-printed bone, with mechanical properties depending on topology and size effect [37,38,39]. Moreover, this kind of predictive tool is important to forecast the expected mechanical behaviour of the scaffold on-site, using an optimizable geometry or bulk material.
In order to do it with FEM, there are two basic steps required. The first one is to represent every scaffold as an equivalent finite element, characterized by an approximately cubic shape with eight nodes, each being able to undergo a displacement in three degrees of freedom. The second step is to discretize the bone geometry in a compatible mesh, made of such solid finite elements, whose material properties are defined to match the mechanical behaviour of the scaffold (Figure 2A).
Accordingly, a flowchart of the required works is presented in Figure 2B. It starts with the scaffold design at iteration i and the material tests for mechanical characterization under tension and compression. Then, an FEM of the bone scaffold is completed, and a monotonic compression test is simulated in that model. Thus, the results from the simulation can be compared to an actual scaffold compression test to calibrate and validate the FE model. This last step is required because the FEM is based on a theoretical or idealized geometry, whereas the real scaffolds manufactured at such scales present faults and uncertainties that need to be taken into account.
For the characterization of additive manufactured scaffolds, an FEM and support methodology has been carried out in this study, whose main objective is to be able to transfer the global mechanical properties of the scaffolds to solid eight-node finite elements, which will be used in FEM analyses. This will enable mechanical calculations to be addressed at a macroscopic level (complete bones, etc.), minimizing the computational costs, because it will not be necessary to take into account, in these models, the exact configuration of the scaffolds, but their equivalent mechanical properties. Additionally, this FEM of additive manufactured scaffolds has different properties, as equivalent finite elements can be further used to model mechanical property gradients simply by progressively changing such finite element material properties accordingly along a certain direction.

2. Materials and Methods

2.1. Materials

The bulk materials composing the scaffolds for bone tissue engineering were selected from amongst recently investigated new formulations based on combinations of a well characterized synthetic block copolymer of poly(ethylene oxide terephthalate)/poly(butylene terephthalate) (PEOT/PBT) with a diverse filler concentration [32,40,41]. For the improvement of the mechanical properties, the scaffolds produced from PEOT/PBT, including 45% of nano hydroxyapatite (nanoHA), were chosen due to their superior properties (Figure 3) and, therefore, were also used for the FE modelling. The mean mechanical behaviour obtained from the stress–strain curves was used (Figure 4). The stress–strain curves were obtained from the mechanical testing of the bulk materials according to the ASTM F2027-08 standard guide [42], prescribing ASTM D638-08 and ASTM D695-10 for tensile and compressive properties, respectively [43,44], and according also to international standard ISO 527 [45,46,47] (Figure 4). A Poisson’s ratio of 0.4 was taken as usual practice for polymers, although there is some evidence that it could be higher in the case of PEOT/PBT [16]. The effect of the bulk material’s Poisson’s ratio on scaffold mechanical properties was studied in a sensitivity analysis.
For the scaffolds, monotonic compression tests, according to ASTM F2150-02 and ISO 604, were carried out on cylindrical scaffolds [48,49] to study the effect of several parameters, namely strand distance, the diameter of the extrusion needle (which determines the scaffold fibre diameter), layer height and the structure of the scaffold (external geometry, fibre position, scaffold diameter, scaffold height). Figure 5).
Three variations of the scaffold were produced and tested, with three specimens each. The variations tested involved a decrease in the strand distance, meaning a denser, less porous scaffold, with a smaller inner volume of void spaces. The geometry of the different scaffold variations is summarized in Table 1.
The monotonic compression tests of the scaffolds showed a good correspondence between samples of the same type within the elastic range, with a clear improvement in the mechanical properties as scaffold porosity decreases and scaffold density is, thus, increased (Figure 6). All scaffolds showed linear elastic behaviour under 10% strain and good reproducibility of stress–strain behaviour, supporting their suitability for load bearing bone tissue engineering, i.e., the intended use. Nevertheless, the differences in the strain–strain relationship of samples of the same type, evidenced in the curves, which were much more pronounced in the 3D-printed scaffolds than in the equivalent testing of bulk materials (Figure 4), indicate subsequent result scattering due to printing accuracy issues. Therefore, the only way to improve the result regularity on the scaffolds is to improve the manufacturing accuracy and quality control, where possible, whose only alternative is to deal with uncertainty, taking it into account during the design and calculation stage for a patient-customized scaffold.
The most relevant values of the elastic range are the Young’s modulus, defined as the slope of the initial linear relationship between the stress and strain, and the yield strength, defined as the stress where the linear range ends and the non-linear behaviour continues developing plasticity. Such values are summarized in Table 2, with corresponding Mean and Standard Deviation (SD) values from the tests of three replicates each.
As a final comment, since this manuscript is about FEM of bone scaffolds, the results of monotonic compression tests on bulk materials and scaffolds are, respectively, just an input to set the material properties in a finite element model and a comparison basis for the output for validation. Thus, since these are boundary conditions of the subject matter of the manuscript, required to reproduce FEM results, they are presented in this section.

2.2. Methods

Several numerical models have been developed in order to simulate the behaviour of the three cylindrical scaffolds mentioned in Section 3.1 when subjected to a monotonic compression test. To this objective, ANSYS v17 finite elements software has been employed.
The numerical models developed within the frame of this study follow the geometrical characteristics defined in Table 2 in terms of fibre diameter (0.25 mm), strand distance (1.25/1.00/0.75 mm) and layer height (0.20 mm). A diameter of 4 mm has also been considered for the samples. Regarding the total height of the scaffolds, due to geometrical considerations to enable the idealized interlayer link modelling, there is little difference between the real height (4.00 mm) and the height considered in the numerical models (4.05 mm).
Considering the double symmetry, not only of the geometry of the scaffolds, but also of the loadings to be considered, just one quarter of the scaffolds has been represented in the numerical models, as shown in Figure 7.
For the modelling of the fibres that conform to the scaffold, only one type of element (SOLID187) has been used. The SOLID187 element is a higher order 3D element with a quadratic displacement behaviour that is well suited to modelling irregular meshes. This element is defined by 10 nodes (tetrahedral element with nodes at vertices and mid edges) with three degrees of freedom at each node (translations in the nodal X, Y and Z directions). Moreover, for the modelling of the contacts, two additional types of elements have been automatically introduced by ANSYS, i.e., CONTA174 and TARGE170. CONTA174 is used to represent contact and sliding between 3D “target” surfaces and a deformable surface defined by this element, whereas the TARGE170 element is used to represent various 3D “target” surfaces for the associated CONTA174 contact elements. Hence, for demonstrative purposes, a detail of the meshing for G25_nanoHA_1.25 FEM with such finite elements is shown in Figure 8.
Regarding the boundary conditions, as the fibre diameter (0.25 mm) is larger than the layer height (0.20 mm), there will be stacking between perpendicular fibres. In these numerical models, it has been considered that there is a perfect bonding between fibres, as shown in Figure 9a. Additionally, a lower horizontal plane has been included in these numerical models to limit, by means of a contact, the vertical displacement of the lower fibres. This plane presenting infinite stiffness and its displacements are totally constrained, see Figure 9b. Finally, an upper horizontal plane has also been included in these numerical models in order to impose, by means of a contact, vertical displacement in the upper fibres. This plane presents infinite stiffness too and its lateral displacements are totally constrained (Figure 9c).
Finally, regarding the actions and type of analysis, in order to simulate a monotonic compression test, a vertical displacement has been imposed on the upper plane of the numerical models. This vertical displacement will be transferred to the scaffold by means of the contact defined between the upper plane and the upper fibres. Due to the negligible influence on the results, the self-weight of the scaffolds has not been considered in this study. Additionally, the need for evaluating the non-linearities related to both the material and the existence of contacts, has demanded the resolution of the numerical models developed in this study by means of a non-linear static analysis. In this resolution, the hypothesis of small displacements has also been taken into account.

3. Results

3.1. G25_nanoHA 1.25

In Figure 10a, the stress–strain curve for the G25_nanoHA_1.25 scaffold is shown. The strain (ε = Δl/l) is then calculated as the displacement of the upper plane divided by the total height of the scaffold, 4.05 mm. Additionally, the stress is calculated as the force to be applied in order to obtain a certain displacement of the upper plane divided by the cross-sectional area of the scaffold (π·φ2/4 = 4π mm2), where φ is the diameter of the cylindrical scaffold. In this regard, this stress can be considered as an “apparent stress”. This is a very important difference, because the stress–strain relationship derived by this way will differ from the actual stresses obtained using the FE model, depicting stresses of every prismatic solid FE according to a colour scale.
Thus, defining the elastic modulus of the scaffold E as the stress divided by the strain for a strain value of 0.03 (within this range, the stress–strain relationship is almost linear) leads to the following Equation (1):
E = 0.1955 0.03 = 6.60   MPa
Additionally, in the following Figure 10b, the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown. Von Mises stress has been chosen for comparison because it is able to combine stresses in several axes in a single stress for comparison. Despite the fact there is a clearly dominant stress direction, since this is a simulated monotonic uniaxial compression test, it is combined with tangential and normal stresses in other directions at the spots where the fibres are crossing each other; therefore, simply plotting the principal stresses will underestimate the real triaxial stress state. As reflected in that figure, loads are mainly transferred through the columns formed by the intersections between perpendicular fibres.
It is noteworthy to remark the differences between the stresses of a stress–strain relationship (Figure 10a) and those depicted in the FEM (Figure 10b). The first going up to 0.5 MPa, while the second is going up to 13.2 MPa at substep 150. This is because, while the first is the quotient between the applied force divided by the cylinder area, the second is the actual stress at every SOLID 187 FE.

3.2. G25_nanoHA_1.00

Analogously, the same procedure is applied to the G25_nanoHA_1.00 numerical model. Accordingly, its stress–strain curve is shown in the following Figure 11a, with the corresponding Von Mises equivalent stress in Figure 11b. Moreover, the elastic modulus of the scaffold E, as the stress divided by the strain for a strain value of 0.03 leads to the following Equation (2):
E = 0.2342 0.03 = 7.90   MPa
It is noteworthy to remark the differences between the stresses in a stress–strain relationship (Figure 11a) and those depicted in the FEM (Figure 11b). The first going up to 0.28 MPa, while the second is going up to 20.4 MPa at substep 150. This is because, while the first is the quotient between the applied force divided by the cylinder area, the second is the actual stress at every SOLID 187 FE.

3.3. G25_nanoHA_0.75

Finally, the stress–strain curve of the G25_nanoHA_0.75 is shown in Figure 12a, while the corresponding equivalent Von Mises Stress is depicted in Figure 12b. The Young’s modulus E is then calculated as the slope of the stress–strain curve at the point with 0.03 strain, as derived in the following Equation (3):
E = 0.2342 0.03 = 15.49   MPa
It is noteworthy to remark the differences between the stresses in the stress–strain relationship (Figure 12a) and those depicted in the FEM (Figure 12b). The first going up to 0.28 MPa, while the second is going up to 13 MPa at substep 150. This is because, while the first is the quotient between the applied force divided by the cylinder area, the second is the actual stress at every SOLID 187 FE.

4. Discussion

The results obtained using an FEM of the different scaffolds G25 nanoHA 1.25, 1.00 and 0.75 are compared to the corresponding monotonic compression tests on the Young’s modulus basis (Table 3). As can be seen in all the cases, the finite element model with idealized geometry showed a Young’s modulus higher than the one measured using direct testing. Moreover, this improvement was higher than the mean value plus 2–3 times the standard deviation of each scaffold. In the case of G25 nanoHA 1.25, the relationship of the FEM value to the mean measured value was 6.60/2.25 = 2.93 times higher, while in the case of the other two, G25 nanoHA 1.00 and G25 nanoHA 0.75, it was 1.54 and 2.45, respectively.
Therefore, after repeating the FEM analysis several times, changing mesh sizes, contact definition and calculation steps and looking at the results, the difference was so high that it could not be explained by result scattering only. In fact, there are several assumptions on the FEM analysis on an idealized scaffold model that are not realistic at all, each contributing to these differences. For instance, a nominal fibre diameter, fibre position and straightness or strand distance are perfect in the model, but in the real world, with a 3D-printing process at such minimal scales, these parameters are impossible to be guaranteed and difficult to be controlled in real time. Thus, such imperfections are sources of uncertainty in terms of mechanical behaviour and the only way to deal with this is to set some execution tolerances and apply a design factor to the idealized geometry, as is common practice in structural engineering to deal with such execution imperfections.
Hence, the first step is to identify the main prima facie uncertainty sources. As such unintentional imperfections are impossible to reproduce, isolated from other errors and with accuracy in real bone scaffolds for testing, the best approach is to model and simulate such imperfections and see their influence on the results. In order to be able to simulate it, each uncertainty source is considered as an independent random variable, whose result is expected to be within a certain interval, with a reasonable grade of reliability, considering an according execution control procedure to guarantee some defined tolerances. Finally, all the variables need to be combined in a single design safety factor considering the global uncertainty, since each one will contribute to it to a certain degree.
As a last comment, a bone scaffold presents time-dependent properties in vivo. In fact, it undergoes deterioration processes related to fatigue, because a bone suffers from cyclic loadings and equivalent biological remodelling acting such as corrosion, since the scaffold material is progressively removed and substituted by bone tissue and this points for an initial deterioration up to a valley point where it starts to recover mechanical properties. Accordingly, the best approach is to develop a simplified model to take into account these two deterioration processes during service life, such as those presented by Calderon Uriszar-Aldaca et al. [50,51,52], developing an additional correction factor to take these into account.

4.1. Sensitivity Analysis

For the sensitivity analysis, the G25_nanoHA_1.25 numerical model under monotonic compression is taken as the base reference model, since it showed the highest difference between the tests and the idealized model, which is already described and whose results are presented in Section 3.1. Hence, the base characteristics and corresponding variations are presented in Table 4.

4.1.1. Poisson’s Ratio of the Base Material

As mentioned in Section 3.1 materials, a Poisson’s coefficient of 0.4 has been considered for the base material as common practice for polymeric materials. Nevertheless, as no specific test has been carried out in the frame of the project in order to determine this value accurately, it can be considered as a source of uncertainty.
Accordingly, three simulations have been carried varying the Poisson’s ratio with 0.4, 0.375 and 0.35 as inputs for G25_nanoHA_1.25 scaffold finite element model. Then, the stress–strain relationship and the elastic modulus has been obtained following the same procedure to see how this variation influences the results.
Thus, Figure 13 shows the stress–strain relationship when varying the Poisson’s ratio and the elastic modulus obtained as the slope of such stress–strain relationship, the curves are so close that they stack over each other. Moreover, Table 5 summarize the results. According to these results, the influence is almost negligible, i.e., after increasing the Poisson’s ratio of the base material, only a slight increase in the elastic modulus of the scaffold can be appreciated, lowering only up to 99.09% in stiffness.

4.1.2. Position of the Fibres within the Scaffold

According to Section 3 on results, the compression loads are mainly transferred to the basement through the columns formed by the intersections between the perpendicular fibres of the scaffold (Figure 10). Hence, due to the relatively low value of the sample diameter in comparison with the strand distance, the relative position of the fibres within the scaffold will affect the number and location of these columns, having some influence on the compressive stiffness of the scaffold, which can happen as a result of punching position variability.
As shown in Figure 14, two different configurations have been analysed. In the first one, at each layer of the scaffold, there is a fibre located in the relevant plane of symmetry. On the other hand, in the second configuration, at each layer of the scaffold, the relevant plane of symmetry is in the centre of the gap between two fibres.
Therefore, a detail of the geometry of the numerical models used for both configurations is shown in the following Figure 15:
The influence of the relative position of the fibres within the scaffold on the stress–strain relationship of the G25_nanoHA 1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 16).
The values of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold of both configurations are shown in Table 6 (the compressive elastic modulus of the scaffold has been defined as the stress divided by the strain for a strain value of 0.03).
According to these results, there is some influence of the relative position of the fibres within the scaffold on the compressive elastic modulus of a G25_nanoHA_1.25 scaffold, with configuration one being stiffer than configuration two, only reaching 94.85% of the stiffness. Finally, the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 17 for both the configurations.

4.1.3. Fibre Diameter

As already disclosed in Table 1, scaffolds have been printed using a theoretical fibre diameter of 0.25 mm. Nevertheless, the accuracy of the printer in terms of deposition speed, positioning, temperature and material viscosity and variations in the material flow rate can lead to the real value of the fibre diameter being different to the theoretical one.
The influence of the fibre diameter on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed by means of FEM (Figure 18).
The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the fibre diameter is shown in Table 7.
The cross-sectional area of the columns formed by the intersections between perpendicular fibres of the scaffold will depend on the fibre diameter (the bigger the fibre diameter, the bigger the cross-sectional area of these columns). As the compression loads applied on the scaffold will be mainly transferred through these columns, the fibre diameter will have a great influence on the compressive elastic modulus of the scaffold, as can be seen in Table 7 and Figure 18.
For two different values of the fibre diameter (φfibre = 0.235 mm and φfibre = 0.220 mm), the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 19, for φfibre = 0.250 mm (see also Figure 10).

4.1.4. Strand Distance

The G25_nanoHA_1.25 scaffold samples have been printed using a theoretical strand distance of 1.25 mm. Nevertheless, the accuracy of the positioning system of the printer can lead to a real value for the strand distance different from the theoretical one.
The influence of the strand distance on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM and the results of this analysis are shown Figure 20.
The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the strand distance is shown in Table 8.
According to these results, when the variation of the strand distance is low, their influence on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load is almost negligible, i.e., increasing the strand distance causes only a slight decrease in the elastic modulus of the scaffold, up to 99.55% of the original stiffness.
As the compression loads applied on the scaffold will be mainly transferred through the columns formed by the intersections between the perpendicular fibres of the scaffold, variations in the strand distance will have a significant impact on the compressive elastic modulus of the scaffold, only when these variations lead to an increase/decrease in the number of those intersections.
For two different values of the strand distance (d = 1.20 mm and d = 1.30 mm), the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 21, for d = 1.25 mm (see also Figure 10).

4.1.5. Layer Height

The tested G25_nanoHA_1.25 scaffolds have been printed using a theoretical layer height of 0.20 mm. The influence of the layer height on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 22).
The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the strand distance is shown in Table 9.
The cross-sectional area of the columns formed by the intersections between the perpendicular fibres of the scaffold will depend on the layer height; the bigger the layer height is, the lower the cross-sectional area of these columns is. As the compression loads applied on the scaffold will be mainly transferred through these equivalent columns, the layer height will have a great influence on the compressive elastic modulus of the scaffold, as can be seen in Table 9 and Figure 22. For two different values of the layer height (h = 0.215 mm and h = 0.230 mm), the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 23, for h = 0.200 mm (see Figure 10).

4.1.6. Sample Diameter

In the frame of this study, monotonic compression tests and dynamic compression tests have been carried out on cylindrical samples with a theoretical diameter of 4 mm, extracted from G25_nanoHA_1.25 scaffolds. Nevertheless, the accuracy of the tools used to extract the samples from the printed scaffolds can lead to a real value for the sample diameter different from the theoretical one. The influence of the sample diameter on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 24).
The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the sample diameter is shown in Table 10.
According to these results, when the variation of the sample diameter is such that the number of columns formed by the intersections between perpendicular fibres of the scaffold is not affected, their influence on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load is negligible. On the other hand, if the actual sample diameter is used to determine the stress on the sample, a certain influence of this parameter on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load can be observed, as can be seen in Table 10 and Figure 24.
These results highlight that, due to the relative low value of the sample diameter in comparison with the strand distance, the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load is very sensitive to the diameter of the tested sample. Nevertheless, this influence is expected to decrease as the sample diameter increases.

4.1.7. Deflection of Fibres during the Printing Process

Due to the low mechanical properties of the fused material, the fibres tend to bend and lose straightness during the printing process at middle span, between rigid or supporting nodes, such as the intersections between the perpendicular fibres, as can be appreciated in Figure 25.
The influence of the deflection of the printed fibres on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (in Figure 26).
The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the deflection of the printed fibres is shown in Table 11.
According to these results, the influence of the deflection of the printed fibres on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load is very low, i.e., increasing the deflection of the printed fibres means a slight decrease in the elastic modulus of the scaffold up to 97.12% of the original stiffness. It should be noted that, in this case, the compression load is parallel to the columns formed by the intersections between the perpendicular fibres. For compression loads parallel to any of the fibre directions, the influence of the deflection of the printed fibres on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load is expected to be much higher.
For two different values of the deflection of the printed fibres, f = 0.06 mm and f = 0.12 mm, the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 27; for no deflection, see Figure 10.

4.1.8. Straightness of the Columns Responsible for Load Transmission

G25_nanoHA_1.25 scaffolds have been printed using a theoretical strand distance of 1.25 mm. However, the accuracy of the positioning system of the printer can lead to a real value for the strand distance different from the theoretical one.
Hence, the influence of the strand distance on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load was analysed in Section 4.1.4 on strand distance. In that case, it was supposed that there was some variation in the strand distance, but it was also considered that this variation remained constant within the whole scaffold. Under these circumstances, the columns formed by the intersections between the perpendicular fibres, which are responsible for the transmission of the compression load applied on the scaffold, remained vertically straight. Nevertheless, the x and y positioning in any printer present deviations; therefore, the real positions are x ± Δx and y ± Δy, meaning that some executions could present alternate positioning or show a displacement tendency.
In this section, on the other hand, the influence of the strand distance on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has also been analysed. In this case, this variation in the strand distance affected the straightness of those columns.
Thus, the following two situations have been analysed:
  • Alternate variations in the strand distance (situation one);
  • Cumulative variations in the strand distance (situation two).
In the first situation, there is a certain eccentricity e between the theoretical location of the fibres and their actual location. In addition, it has been considered that the fibre corresponding to a certain layer is located on the opposite side, regarding the theoretical location of the fibres, of that of the fibres corresponding to the preceding and the following layers, see Figure 28.
A detail of the geometry of the numerical model used to simulate this situation is shown in Figure 29. In this situation, the following two values have been considered for the deviation of the fibres: e = 0.025 mm and e = 0.05 mm.
In the second situation, there is also a certain eccentricity between the theoretical location of the fibres and their actual location. On the other hand, in this case, the deviation is cumulative, as shown in Figure 30.
A detail of the geometry of the numerical model used to simulate this situation is shown in Figure 31. In this situation, two values have also been considered for the deviation of the fibres, e = 0.005 mm and e = 0.009 mm.
In both situations, it has been considered that the deviation between the theoretical location of the fibres and their actual location affects the fibres oriented in both directions.
As can be seen in Figure 29 and Figure 31, in these cases, there is no any symmetry in the geometry of the scaffolds; therefore, the whole scaffold has been represented in the numerical models. Therefore, in order to avoid lateral rigid body motions, which may impede the convergence of the numerical models, a friction coefficient of 0.3 has been considered in the contact between the upper plane and the upper fibres and the contact between the lower plane and the lower fibres.
The influence of the deviation between the theoretical location of the fibres and their actual location, affecting the straightness of the columns formed by the intersections between the perpendicular fibres, on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 32).
Accordingly, the values of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis, for the deviation between the theoretical location of the fibres and their actual location, are shown in Table 12 for situation one and Table 13 for situation two.
As can be seen in Table 12 and Table 13, the compressive elastic modulus of the scaffold with no deviation (6.70 MPa) differs from the compressive elastic modulus of the reference scaffold considered in the previous sections (6.60 MPa). This difference lies in the fact that, as mentioned, the numerical models used in both cases differ regarding the friction coefficient used for the contacts between the upper and lower planes and the scaffold fibres.
According to these results, when errors are cumulative (situation two), there is no significative influence, for the range of values considered, of the deviation between the theoretical location of the fibres and their actual location on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load. However, when errors alternate with respect to the theoretical location of the fibres (situation one), the influence of the deviation between the theoretical location of the fibres and their actual location on the stress–strain relationship of a G25_nanoHA 1,25 scaffold subjected to a compression load may be considerable, because the straightness of the columns responsible for load transmission is much more affected in such situation.

4.1.9. Existence of Broken Fibres

The G25_nanoHA_1.25 scaffold samples for testing have been printed using a theoretical fibre diameter of 0.25 mm, leading to relatively fragile fibres that can suffer damage during the printing process, the cooling and polymerization process or the tooling of the samples.
Thus, the influence of the existence of broken fibres on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM, having considered the following alternatives:
  • No broken fibres within the scaffold—Situation one.
  • One broken fibre within the scaffold (breakage located in one intersection between perpendicular fibres)—Situation two.
  • Two broken fibres within the scaffold (breakages located in two different intersections between perpendicular fibres)—Situation three.
  • One broken fibre within the scaffold (breakage located far from any intersection between perpendicular fibres)—Situation four.
  • Two broken fibres within the scaffold (breakages located far from any intersection between perpendicular fibres)—Situation five.
All the fibre breakages considered in this study had a length of 0.15 mm and their location was randomly selected.
As in the previous section, in this case there is no symmetry in the geometry of the scaffolds; therefore, the whole scaffold has been represented in the numerical models. Therefore, in order to avoid lateral rigid body motions, which may impede the convergence of the numerical models, a friction coefficient of 0.3 has been considered in the contact between the upper plane and the upper fibres and the contact between the lower plane and the lower fibres. The results of this analysis are shown in Figure 33.
The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the five different situations considered in this analysis concerning the existence of broken fibres is shown in Table 14.
As can be seen in Table 14, the compressive elastic modulus of the scaffold with no broken fibres (6.70 MPa) differs from the compressive elastic modulus of the reference scaffold considered in the previous sections (6.60 MPa). This difference lies in the fact that, as mentioned, the numerical models used in both cases differ regarding the friction coefficient used for the contacts between the upper and lower planes and the scaffold fibres.
According to these results, when the breakage of the fibres occurs far from any intersection between perpendicular fibres, there is no significant influence on the stress–strain relationship of a G25_nanoHA 1.25 scaffold subjected to a compression load, as the paths for the load transmission (columns formed by the intersections between perpendicular fibres) are not affected.
On the other hand, some influence on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load can be seen when the breakage of the fibres is located in the intersection between perpendicular fibres. In this case, it is expected that the longer the breakage is and the higher the number of affected fibres is, the bigger the influence is.

4.2. Design Safety Factor to Consider Uncertainty

In view of the Section 4.1 outcome, it is clear that the great differences in terms of mechanical behaviour summarized in Table 3 could be explained by the printing process uncertainty. Any model will be performed according to an idealized geometry, conditions and material properties but, in the end, assuming this idealized behaviour is unsafe because 3D printing bone scaffolds at such small-scale leads to several manufacturing mistakes, which become uncertainty sources.
Nevertheless, although it is not possible to forecast exact predictions on mechanical behaviour while unable to perform a perfect 3D-printed bone scaffold, it is possible to perform enough safe ones. This procedure is analogous to that used in purely structural or mechanical engineering for other non-medical applications, where the forecasts are performed running idealized models with reduction factors to consider execution mistakes, material properties uncertainties, lack of homogeneity, etc.
Thus, such reduction factors need to be balanced to be conservative enough but not too penalizing for the structure and the building process. Therefore, several tolerances are defined for each execution variable that need to be checked and controlled during structure execution or, analogously in this case, for bone scaffold 3D printing. Hence, each variable is defined by a mean-measured value, μ, and a standard deviation, σ, assuming a Gaussian distribution according to central limit theorem, since each variable is dependent on several uncontrolled factors, and defining the check and control intensity in a way that there is a small probability, typically 5%, of being under a lower boundary limit or exceeding an upper bound with a reliability level typically of 75%. As an additional comment, these fifth percentile values are known as characteristic values.
Accordingly, each uncertainty source is defined as a random variable x1, x2, …, x9 and, since the influence of each isolated random variable on the elastic modulus is linear dependent and biunivocal. the value of each random variable corresponds with an uncertainty factor, Fi, that is defined as the quotient between the actual elastic modulus and the theoretical one if 3D printed perfectly. For instance, the fibre diameter is the variable x3, the nozzle diameter is 0.25 mm, and it theoretically prints at 0.25 mm but, depending on the printing speed, the fibre diameter could be lower. Thus, if the 3D printing process was perfect, then the fibre diameter would be 0.25 mm and the elastic modulus would be E = 6.6 MPa with a corresponding uncertainty factor F3 = 6.6/6.6 = 1.0. However, if the printing speed was higher, then the fibre diameter would be lower instead and, in the case that it was even 0.22 mm, then the elastic modulus would be E = 3.39 MPa and the corresponding factor F3 = 3.39/6.6 = 0.5136. Thus, because of the linear dependency, F3 is a random variable with a Gaussian distribution matching x3. Then, if we define a diameter tolerance and perform an execution control to guarantee that the fibre diameter is kept between 0.22 and 0.25 mm, with such an intensity that only 5% of the bone scaffolds present fibres below 0.22 mm and only 5% over 0.25 mm. Analogously, a mean, μi, and a standard deviation, σi, can be defined for each isolated uncertainty factor Fi, considering that in any Gaussian distribution, 90% of the probability is within the mean minus 1.6449 times the standard deviation and plus it. Hence, Table 15 summarizes the mean, μi, and standard deviation, σi, of each isolated uncertainty factor Fi, with corresponding limit values, Fi,min and Fi,max, for the interval.
Thus, the next Figure 34 shows the probability density functions fi(xi) of each Gaussian random variable thus defined. For the sake of simplicity, each variable xi directly represents the uncertainty factor Fi = E(xi)/E(xtheoretical) defined as the quotient between the actual elastic modulus and the theoretical one if 3D printed perfectly.
The mathematical expression of each Gaussian probability density function corresponds to Equation (4), as follows:
f i ( x i ) = 1 σ i · 2 π · e ( x i μ i ) 2 2 · σ i 2
Hence, considering each isolated uncertainty variable as independent from each other, the probability density function of the multiplication is the multiplication of the independent probability density functions, see Equation (5) as follows:
f ( x 1 · x 2     x 9 ) = f 1 ( x 1 ) · f 2 ( x 2 )     f 9 ( x 9 )
Therefore, to derive the probability density function of the multiplication it is required a variable change. Accordingly, such variable change is shown in Equation (6), where u is the multiplication of each N = 9 variables xi and there are other N-1 variables zi.
{ u = x 1 · x 2     x 9 z 1 = x 1 z 2 = x 2 z 8 = x 8
Accordingly, the Jacobian determinant to execute the variable change is derived as the partial derivation of the variable multiplication by each variable, see Equation (7).
( x 1 · x 2     x 9 ) ( u , z 1 · z 2     z 8 ) = | x 1 u x 1 z 1 x 1 z 8 x 2 u x 2 z 1 x 2 z 8 x 9 u x 9 z 1 x 9 z 8 | = | 0 1 x 1 z 8 0 0 x 2 z 8 1 z 1 · z 2     z 8 u z 1 2 · z 2     z 8 u z 1 · z 2     z 8 2 | = 1 z 1 · z 2     z 8
Thus, the probability density function of the nine variables, once changed, is shown in Equation (9), as follows:
( u , z 1 , , z 8 ) = f 1 ( z 1 ) · f 2 ( z 2 )     f 9 ( u z 1 · z 2     z 9 ) · 1 z 1 · z 2     z 8
In order to obtain the probability density function, the marginal distribution for variable u, which is obtained by the direct integration of the remaining variables z1 to z8, is required. See Equation (9) as follows:
( u ) =   R 8 [ f 1 ( z 1 ) · f 2 ( z 2 )     f 9 ( u z 1 · z 2     z 9 ) · 1 z 1 · z 2     z 8 ] · dz 1 · dz 2     dz 8
As a last comment on Equation (9), it is worthy to remember that each fi(xi) is a normal Gaussian distribution according to Equation (4) and that, despite being multiplied by the Jacobian fraction at the end, it is not directly integrable. Therefore, any solution can be obtained only by complex numerical means. Nevertheless, there is another simpler, more practical approach, considering that the distribution of the multiplication of Gaussian distributions is a Gaussian distribution itself. Then, it is possible to derive the mean and the standard deviation simply by taking into account their properties and their relationship with the expected value, E(x), and variance, Var(x).

4.2.1. Derivation of the Global Mean

The expected value corresponds with the first order moment of the probability distribution. Therefore, the expected value of the multiplication of variables is the ninth integration of the probability density function of the multiplication multiplied by each variable xi, between −∞ and ∞. See Equation (10), which is as follows:
E [ x ^ 1 · x ^ 2     x ^ 9 ] =   R 9 [ x 1 · x 2     x 9 · f x ^ 1 · x ^ 2     x ^ 9 ( x 1 · x 2     x 9 ) · dx 1 · dx 1     dx 9 ]
Nevertheless, as the probability density function of the multiplication for a set of independent variables is the product of the isolated probability density function of each variable, see Equation (5), then the Equation (10) can be developed into Equation (11), which is as follows:
E [ x ^ 1 · x ^ 2 x ^ 9 ] =   R 9 [ x 1 · x 2 x 9 · f x ^ 1 ( x 1 ) · f x ^ 2 ( x 2 ) f x ^ 9 ( x 9 ) · dx 1 · dx 1 dx 9 ]
Accordingly, as every integral is independent from each other, then they can be regrouped as in the following Equation (12), presenting them as the product of integrals:
E [ x ^ 1 · x ^ 2 x ^ 9 ] = x 1 · f x ^ 1 ( x 1 ) · dx 1 · x 2 · f x ^ 2 ( x 2 ) · dx 2     x 9 · f x ^ 9 ( x 9 ) · dx 9
Thus, remembering that the mean is precisely the first order moment of a probability distribution, see Equation (13), which is as follows:
x i · f x ^ i ( x i ) · dx i = μ i
Then, by direct substitution, the mean of any variable multiplication is the product of the isolated variable means, see Equation (14), which is as follows:
E [ x ^ 1 · x ^ 2 x ^ 9 ] = μ 1 · μ 2 · μ 3     μ 9

4.2.2. Derivation of the Global Standard Deviation

Regarding the variance of the variable multiplication, it corresponds to the expected value of the squared product minus the squared expected value of the product, see (15), which is as follows:
VAR ( x 1 · x 2 · x 3     x 9 ) = E [ ( x 1 · x 2 · x 3 x 9 ) 2 ] E ( [ x 1 · x 2 · x 3     x 9 ] ) 2
Now, making use of the expected value properties, the expected value of the product of variables is the multiplication of the expected value of such variables. Therefore, substituting in Equation (15) there it comes Equation (16), which is as follows:
VAR ( x 1 · x 2 · x 3     x 9 ) = E [ x 1 2 ] · E [ x 2 2 ]     E [ x 9 2 ] ( E [ x 1 ] ) 2 · ( E [ x 2 ] ) 2     ( E [ x 9 ] ) 2
Now, considering that the expected value of a variable is the mean value E[xi] = μi, the expected value of the squared variable corresponds to Equation (17), which is as follows:
E [ x i 2 ] = ( μ i 2 + σ i 2 )
Then, by a simple substitution in Equation (16), Equation (18) is then derived as follows:
VAR ( x 1 · x 2 · x 3     x 9 ) = ( μ 1 2 + σ 1 2 ) · ( μ 2 2 + σ 2 2 )     ( μ 9 2 + σ 9 2 ) μ 1 2 · μ 2 2     μ 9 2

4.2.3. Derivation of the Safety Factor

Thus, according to previous Equations (14) and (18), when substituting the mean and standard deviation values for the uncertainty factors, taken as isolated variables, which are summarized Table 15, the mean and standard deviation of the product of variables is derived. See Equations (19) and (20), which are as follows:
μ = E [ x ^ 1 · x ^ 2 x ^ 9 ] 0.4578
σ = VAR ( x 1 · x 2 · x 3     x 9 ) 0.1346
Therefore, with an execution control of each isolated uncertainty source, keeping 95% of the bone scaffolds within the required intervals for the fibre diameter, layer height, etc., then the combined outcome in terms of the global uncertainty factor will show a Gaussian distribution (Figure 35). In this figure, the mean value is 0.4578, while the characteristic value corresponding to fifth percentile (shadowed area) is 0.2364.
Consequently, there are two reduction factors derived from uncertainty to be applied to the elastic modulus directly obtained from the bulk material as correction factors for the scaffolds modelling. The first is derived from the mean value and shall be used for FEM simulations where the accuracy in terms of deformations and displacements are the main outcome, and the second is derived from the characteristic value for conservative forecasts that shall be used for FEM simulations where the main outcome is a safe enough prediction in terms of stress or an upper bound of deformations. Thus, the design value for the global elastic or Young’s modulus after printing Ed can be derived from the elastic modulus obtained from models with idealized or theoretical conditions, EM, and material input obtained from bulk material testing, but with a reduction material printing factor as a design safety factor γMP, see Equation (21). Thus, for accurate simulations regarding deformations it takes 2.2 as value and for safer conservative simulations regarding stresses or the upper bound of local deformations it takes 4.25 instead, see Equations (22) and (23), which are as follows:
E d = E M γ MP
E d = 0.4578 · E M E M 2.2
E d = 0.2364 · E M E M 4.25
Additionally, this design safety factor can be lowered if the execution control becomes more intense; thus 99% of the bone scaffolds presents the variables within the defined tolerance limits instead 95%, or if the intensity remains the same but the tolerance limits are lowered to ensure that geometry deviations are mitigated and derived, the uncertainty factors are lowered accordingly. The fibre diameter, layer height and straightness of the columns are key uncertainty sources affecting the bone scaffold behaviour. Therefore, the best strategy to improve the resulting scaffold behaviour is to improve the execution control and reduce the tolerances in these three variables.

4.2.4. Use Case

For instance, let us have a use case analogous to that presented in Figure 1, where a hollow cylindric scaffold made of G25 nanoHA 1.25, G25 nanoHA 1.00 or G25 nanoHA 0.75 is substituting some damaged length of a femur for tissue regeneration, as representative of a uniaxial loading case. The specific patient body weight is around 75 kg; therefore, the scaffold section must face F = 735 N of loading peaks. If the cylindric scaffold length was 20 mm, the external diameter was 50 mm and the inner diameter was 10 mm, then the cross-sectional area would be A = 1885 mm2. Thus, the axial apparent stress of the scaffold would be σx = F/A = −0.39 MPa, with σy = σz = 0 Mpa, since this is a monoaxial loading and a cylindric loading state, and the axial strain εx is, according to Hooke’s Equation (24) and taking into account the Poisson coefficients μyx and μzx, as follows:
ε x = σ x E x μ yx · σ y E y μ zx · σ z E z = σ x E x
Now, using the apparent young modulus obtained from the idealized models EA,FEM already presented in Table 3 without any correction factor, the mean young modulus, Em,d, does not present any change. Thus, the mean strain εm is obtained by direct application of Equation (24). Then, the displacement δ is obtained by multiplying the strain by the scaffold length. Moreover, the real young modulus, EFEM, taking into account the real stresses of solid 187 elements at the scaffold, multiplied by the strain εm considered as uniform, gives the real foreseeable stress, σFEM. This calculation process for each scaffold type is presented in Table 16.
Now, if the design limits were σmax = 40 MPa and δmax = 1.5 mm, then the selected scaffold type could be G25 nanoHA 0.75 or G25 nanoHA 1.25, as they are the only ones fulfilling both requirements regarding the maximum displacement and stress. Nevertheless, if the scaffold was finally made of such materials, considering that the realistic scaffold will present lower stiffness due to several uncertain sources, then the calculation needs to be readjusted, taking into account the mean expectable behaviour of the scaffold (Table 17). This can be made by applying the reduction factor γMP = 2.2 to derive a mean Young’s modulus for design Em,d and repeating the calculations with this adjusted elasticity. As can be seen, by calculating it this way, the only scaffold type left fulfilling both requirements is G25 nanoHA 0.75.
However, while conducting the calculations this way is a correct procedure to forecast displacements depending on a mean scaffold behaviour, it is not a safe enough procedure to consider the mechanical failure due to stress, since a crack starting at certain spot of the scaffold due to a local weakness could propagate, causing a global structural failure. This is the reason why the calculations on the resistance take a characteristic value instead, i.e., the fifth percentile or the value that is exceeded 19 times of every 20 measurements. This guarantees that there is a certain amount of resistance left to compensate for local weaknesses. Moreover, it can be made by applying the reduction factor γMP = 4.25 instead to derive a characteristic Young’s modulus for design Ek,d and repeating the calculations with this adjusted elasticity. As can be seen, calculating it this way discards the three initial designs, since the three scaffolds present a stress over 40 MPa (Table 18).
The outcome of this evidence is that, if we cannot reduce the loadings, we need to increase the effective cross-sectional area to face them. This can be made by increasing the external diameter or reducing the inner diameter of the scaffold, but most of the time it is not possible in a patient customized geometry; therefore, the only alternatives are as follows:
  • Increasing the execution control, enabling to reduce the factor γMP;
  • Improving the material, pushing forward the stress limits and the global stiffness;
  • Increasing the fibre diameter and/or reducing the strand distance, thus increasing the column area.
Reducing the strand distance from 0.75 to 0.5 leads to an increase in the amount of “resisting columns”, from 21 to 45; this results in an expectable stress of 37.47 MPa, taking into account the previous results in columns made of spatially crossed cylinders of the same geometry. Moreover, increasing the fibre diameter instead, to pass from 60.2 to 40 MPa, will require a cross sectional augmentation by a scale factor equal to the square root of the quotient F = (60.2/40)0.5, giving a 0.307-mm diameter. Nevertheless, the diameter augmentation reduces the fibre cross-sectional curvature, increasing contact in a non-linear basis; therefore, it should be enough to increase the diameter up to 0.3 mm.
Thus, this patient will require a G25 nanoHA 0.75 with a fibre diameter of 0.3 mm and a G25 nanoHA 0.5 with fibre diameter of 0.25 mm, both produced from 55% PEOT/PBT + 45% nanoHA material or, finally, a G25 nanoHA 0.75 with 0.25 mm fibre diameter produced from an improved material.

5. Conclusions

Thus, comparing the results of finite element modelling to the bone scaffold compression tests according to considered methodology, and once the uncertainty sources derived from manufacturing process are analysed during discussion, the following conclusions can be summarized:
A methodology to perform FEM of 3D-printed bone scaffolds of any initial bulk material properties and geometry is defined. Moreover, it is experimentally calibrated and validated by direct testing of the scaffolds. The bone scaffold is then represented as an equivalent finite element to enable the FEM of bone tissue discretized with a compatible mesh, enabling the required characteristics to be studied, depending on bone topology.
The methodology is executed for a 55% PEOT/PBT + 45% nanoHA bulk material and three different scaffold geometries. The results coming from the FEM and the monotonic compression test of the bone scaffolds show significative differences, being that the idealized simulation is always better than the actual bone scaffolds. Considering the tests’ mean values and standard deviations, the differences cannot be explained by result scattering only.
A deep look at the scaffolds show several fabrication imprecisions, meaning a resulting geometry far different from the theoretical one considered for FEM, such as different strand distance, layer height, fibre diameter, Poisson’s ratio, fibre position, sample diameter, straightness of the columns, existence of broken fibres. This is due to the small imprecisions adding up layer-by-layer as large clinically relevant scaffolds (i.e., volume > 1 cm3) are being built. Thus, there are up to nine uncertainty sources detected in scaffolds coming from different variables.
Consequently, a sensitivity analysis is performed for each uncertainty source, deriving the influence or uncertainty factor and reducing the expected behaviour of each isolated source in terms of scaffold stiffness. Then, some tolerances are defined for each variable to be controlled during execution, with such an intensity that every variable presents a 95% probability of being within the selected interval in a Gaussian distribution.
The two main effects that can give strong changes to the elastic modulus are the strand overlapping and the fibre diameter. The first can be accommodated by design planning the shifting of the planes while printing, and the second may be modulating the printing speed, these will contribute to palliate the effect and minimize the corresponding design factor. A global design safety factor is defined in a predictive and conservative scenario, combining each isolated uncertainty factor in a single value taken from the mean value and fifth percentile of the Gaussian distribution of the product of each reducing uncertainty factor.
Thus, the aprioristic FEM of a bone tissue region with the inclusion of bone scaffold volume made of a certain bulk material becomes possible simply by using solid 185 equivalent finite elements to mesh it, assigning as material properties the ones derived from the FEM simulation of the monotonic tests of the scaffold, divided by the global design safety factor. This enables the further topological optimization of the scaffold along the bone tissue and the material selection, depending on the requirements.

Author Contributions

Conceptualization and methodology, I.C.-U.-A.; finite element simulations, S.P.; research on bulk materials and selection, A.S. and S.V.; experimental tests and validation, A.S. and S.V.; formal analysis and investigation, I.C.-U.-A. and A.M; providing resources and scaffolds for testing, R.S. and M.C.-T.; data curation, A.S., S.V. and S.P.; writing—original report, S.P.; writing—manuscript, I.C.-U.-A.; writing—review and editing, I.C.-U.-A., R.S. and L.M.; visualization, A.M. and C.M.; supervision, I.C.-U.-A.; project administration, A.S., A.P. and L.M.; funding acquisition, C.M., A.P. and L.M. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the European Union, represented by the European Commission, grant number 685825-FAST-H2020-NMP-2014-2015/H2020-NMP-PILOTS-2015.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the corresponding authors.

Acknowledgments

Authors would like to acknowledge Miguel Calderon Felices for his help and scientific advice regarding the statistical development of the design safety factor.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Stevens, M.M. Biomaterials for Bone Tissue Engineering. Mater. Today 2008, 11, 18–25. [Google Scholar] [CrossRef]
  2. Fernandez-Yague, M.A.; Abbah, S.A.; McNamara, L.; Zeugolis, D.I.; Pandit, A.; Biggs, M.J. Biomimetic Approaches in Bone Tissue Engineering: Integrating Biological and Physicomechanical Strategies. Adv. Drug Deliv. Rev. 2015, 84, 1–29. [Google Scholar] [CrossRef]
  3. Rezwan, K.; Chen, Q.Z.; Blaker, J.J.; Boccaccini, A.R. Biodegradable and Bioactive Porous Polymer/Inorganic Composite Scaffolds for Bone Tissue Engineering. Biomaterials 2006, 27, 3413–3431. [Google Scholar] [CrossRef]
  4. Jakus, A.E.; Rutz, A.L.; Jordan, S.W.; Kannan, A.; Mitchell, S.M.; Yun, C.; Koube, K.D.; Yoo, S.C.; Whiteley, H.E.; Richter, C.-P.; et al. Hyperelastic “Bone”: A Highly Versatile, Growth Factor–Free, Osteoregenerative, Scalable, and Surgically Friendly Biomaterial. Sci. Transl. Med. 2016, 8, 358ra127. [Google Scholar] [CrossRef] [Green Version]
  5. Karageorgiou, V.; Kaplan, D. Porosity of 3D Biomaterial Scaffolds and Osteogenesis. Biomaterials 2005, 26, 5474–5491. [Google Scholar] [CrossRef]
  6. Hutmacher, D.W. Scaffolds in Tissue Engineering Bone and Cartilage. Biomaterials 2000, 21, 2529–2543. [Google Scholar] [CrossRef]
  7. Lee, W.C.; Lim, C.H.Y.X.; Shi, H.; Tang, L.A.L.; Wang, Y.; Lim, C.T.; Loh, K.P. Origin of Enhanced Stem Cell Growth and Differentiation on Graphene and Graphene Oxide. ACS Nano 2011, 5, 7334–7341. [Google Scholar] [CrossRef] [PubMed]
  8. Chung, C.; Kim, Y.-K.; Shin, D.; Ryoo, S.-R.; Hong, B.H.; Min, D.-H. Biomedical Applications of Graphene and Graphene Oxide. Acc. Chem. Res. 2013, 46, 2211–2224. [Google Scholar] [CrossRef] [PubMed]
  9. Raucci, M.G.; Giugliano, D.; Longo, A.; Zeppetelli, S.; Carotenuto, G.; Ambrosio, L. Comparative Facile Methods for Preparing Graphene Oxide-Hydroxyapatite for Bone Tissue Engineering. J. Tissue Eng. Regen. Med. 2017, 11, 2204–2216. [Google Scholar] [CrossRef]
  10. Badar, M.; Rahim, M.I.; Kieke, M.; Ebel, T.; Rohde, M.; Hauser, H.; Behrens, P.; Mueller, P.P. Controlled Drug Release from Antibiotic-Loaded Layered Double Hydroxide Coatings on Porous Titanium Implants in a Mouse Model. J. Biomed. Mater. Res. A 2015, 103, 2141–2149. [Google Scholar] [CrossRef] [PubMed]
  11. Ai, J.; Chen, Y.; Jing, H.; Gao, X.; Liu, J.; Ma, K.; Suo, J.; Yu, S.; Song, S. Dynamic Release of Antibiotic Drug Gentamicin Sulfate from Novel Zirconium Phosphate Nano-Platelets. Sci. Adv. Mater. 2014, 6, 2603–2610. [Google Scholar] [CrossRef]
  12. Rodriguez-Losada, N.; Wendelbo, R.; Garcia-Fernandez, M.; Pavia, J.; Martinez-Montañez, E.; Lara-Muñoz, J.P.; Arenas, E.; Aguirre-Gomez, J.A. Graphene Derivative Scaffolds Facilitate in Vitro Cell Survival and Maturation of Dopaminergic SN4741 Cells. Acta Physiol. 2014, 212, 69. [Google Scholar] [CrossRef] [Green Version]
  13. Fiorillo, M.; Verre, A.F.; Iliut, M.; Peiris-Pagés, M.; Ozsvari, B.; Gandara, R.; Cappello, A.R.; Sotgia, F.; Vijayaraghavan, A.; Lisanti, M.P. Graphene Oxide Selectively Targets Cancer Stem Cells, across Multiple Tumor Types: Implications for Non-Toxic Cancer Treatment, via “Differentiation-Based Nano-Therapy”. Oncotarget 2015, 6, 3553–3562. [Google Scholar] [CrossRef] [Green Version]
  14. Butscher, A.; Bohner, M.; Roth, C.; Ernstberger, A.; Heuberger, R.; Doebelin, N.; von Rohr, P.R.; Müller, R. Printability of Calcium Phosphate Powders for Three-Dimensional Printing of Tissue Engineering Scaffolds. Acta Biomater. 2012, 8, 373–385. [Google Scholar] [CrossRef]
  15. Butscher, A.; Bohner, M.; Hofmann, S.; Gauckler, L.; Müller, R. Structural and Material Approaches to Bone Tissue Engineering in Powder-Based Three-Dimensional Printing. Acta Biomater. 2011, 7, 907–920. [Google Scholar] [CrossRef] [PubMed]
  16. Moroni, L.; de Wijn, J.R.; van Blitterswijk, C.A. 3D Fiber-Deposited Scaffolds for Tissue Engineering: Influence of Pores Geometry and Architecture on Dynamic Mechanical Properties. Biomaterials 2006, 27, 974–985. [Google Scholar] [CrossRef]
  17. Woodfield, T.B.F.; Malda, J.; de Wijn, J.; Péters, F.; Riesle, J.; van Blitterswijk, C.A. Design of Porous Scaffolds for Cartilage Tissue Engineering Using a Three-Dimensional Fiber-Deposition Technique. Biomaterials 2004, 25, 4149–4161. [Google Scholar] [CrossRef] [PubMed]
  18. Tjong, S.C. Structural and Mechanical Properties of Polymer Nanocomposites. Mater. Sci. Eng. R Rep. 2006, 53, 73–197. [Google Scholar] [CrossRef]
  19. Paul, D.R.; Robeson, L.M. Polymer Nanotechnology: Nanocomposites. Polymer 2008, 49, 3187–3204. [Google Scholar] [CrossRef] [Green Version]
  20. Schieker, M.; Seitz, H.; Drosse, I.; Seitz, S.; Mutschler, W. Biomaterials as Scaffold for Bone Tissue Engineering. Eur. J. Trauma 2006, 32, 114–124. [Google Scholar] [CrossRef]
  21. Pavlidou, S.; Papaspyrides, C.D. A Review on Polymer–Layered Silicate Nanocomposites. Prog. Polym. Sci. 2008, 33, 1119–1198. [Google Scholar] [CrossRef]
  22. Alexandre, M.; Dubois, P. Polymer-Layered Silicate Nanocomposites: Preparation, Properties and Uses of a New Class of Materials. Mater. Sci. Eng. R Rep. 2000, 28, 1–63. [Google Scholar] [CrossRef]
  23. Potts, J.R.; Dreyer, D.R.; Bielawski, C.W.; Ruoff, R.S. Graphene-Based Polymer Nanocomposites. Polymer 2011, 52, 5–25. [Google Scholar] [CrossRef] [Green Version]
  24. Liu, Q.; de Wijn, J.R.; van Blitterswijk, C.A. Composite Biomaterials with Chemical Bonding between Hydroxyapatite Filler Particles and PEG/PBT Copolymer Matrix. J. Biomed. Mater. Res. 1998, 40, 490–497. [Google Scholar] [CrossRef]
  25. Liu, Q.; de Wijn, J.R.; de Groot, K.; van Blitterswijk, C.A. Surface Modification of Nano-Apatite by Grafting Organic Polymer. Biomaterials 1998, 19, 1067–1072. [Google Scholar] [CrossRef]
  26. Liu, Q.; de Wijn, J.R.; van Blitterswijk, C.A. Nano-Apatite/Polymer Composites: Mechanical and Physicochemical Characteristics. Biomaterials 1997, 18, 1263–1270. [Google Scholar] [CrossRef] [Green Version]
  27. Liu, O.; Du Wijn, J.R.; Van Blitterswijk, C.A. Intermolecular complexation between peg/pbt block copolymer and polyelectrolytes polyacrylic acid and maleic acid copolymer. Eur. Polym. J. 1997, 33, 1041–1047. [Google Scholar] [CrossRef]
  28. Liu, Q.; De Wijn, J.R.; Bakker, D.; Van Blitterswijk, C.A. Surface Modification of Hydroxyapatite to Introduce Interfacial Bonding with PolyactiveTM 70/30 in a Biodegradable Composite. J. Mater. Sci. Mater. Med. 1996, 7, 551–557. [Google Scholar] [CrossRef]
  29. Munarin, F.; Petrini, P.; Gentilini, R.; Pillai, R.S.; Dirè, S.; Tanzi, M.C.; Sglavo, V.M. Micro- and Nano-Hydroxyapatite as Active Reinforcement for Soft Biocomposites. Int. J. Biol. Macromol. 2015, 72, 199–209. [Google Scholar] [CrossRef]
  30. Nandakumar, A.; Cruz, C.; Mentink, A.; Birgani, Z.T.; Moroni, L.; van Blitterswijk, C.; Habibovic, P. Monolithic and Assembled Polymer—Ceramic Composites for Bone Regeneration. Acta Biomater. 2013, 9, 5708–5717. [Google Scholar] [CrossRef] [PubMed]
  31. Fouad, H.; Elleithy, R.; Alothman, O.Y. Thermo-Mechanical, Wear and Fracture Behavior of High-Density Polyethylene/Hydroxyapatite Nano Composite for Biomedical Applications: Effect of Accelerated Ageing. J. Mater. Sci. Technol. 2013, 29, 573–581. [Google Scholar] [CrossRef]
  32. Sinha, R.; Cámara-Torres, M.; Scopece, P.; Falzacappa, E.V.; Patelli, A.; Moroni, L.; Mota, C. A Hybrid Additive Manufacturing Platform to Create Bulk and Surface Composition Gradients on Scaffolds for Tissue Regeneration. bioRxiv 2020, 165605. [Google Scholar] [CrossRef]
  33. Olivares, A.L.; Marsal, È.; Planell, J.A.; Lacroix, D. Finite Element Study of Scaffold Architecture Design and Culture Conditions for Tissue Engineering. Biomaterials 2009, 30, 6142–6149. [Google Scholar] [CrossRef]
  34. Miranda, P.; Pajares, A.; Guiberteau, F. Finite Element Modeling as a Tool for Predicting the Fracture Behavior of Robocast Scaffolds. Acta Biomater. 2008, 4, 1715–1724. [Google Scholar] [CrossRef]
  35. Sandino, C.; Planell, J.A.; Lacroix, D. A Finite Element Study of Mechanical Stimuli in Scaffolds for Bone Tissue Engineering. J. Biomech. 2008, 41, 1005–1014. [Google Scholar] [CrossRef] [PubMed]
  36. Milan, J.-L.; Planell, J.A.; Lacroix, D. Computational Modelling of the Mechanical Environment of Osteogenesis within a Polylactic Acid–Calcium Phosphate Glass Scaffold. Biomaterials 2009, 30, 4219–4226. [Google Scholar] [CrossRef] [Green Version]
  37. Almeida, H.A. Design of Tissue Engineering Scaffolds Based on Hyperbolic Surfaces: Structural Numerical Evaluation. Med. Eng. 2014, 8, 1033–1040. [Google Scholar] [CrossRef] [PubMed]
  38. Dinis, J.C.; Morais, T.F.; Amorim, P.H.J.; Ruben, R.B.; Almeida, H.A.; Inforçati, P.N.; Bártolo, P.J.; Silva, J.V.L. Open Source Software for the Automatic Design of Scaffold Structures for Tissue Engineering Applications. Procedia Technol. 2014, 16, 1542–1547. [Google Scholar] [CrossRef]
  39. De Amorim Almeida, H.; da Silva Bártolo, P.J. Virtual Topological Optimisation of Scaffolds for Rapid Prototyping. Med. Eng. Phys. 2010, 32, 775–782. [Google Scholar] [CrossRef]
  40. Cámara-Torres, M.; Duarte, S.; Sinha, R.; Egizabal, A.; Álvarez, N.; Bastianini, M.; Sisani, M.; Scopece, P.; Scatto, M.; Bonetto, A.; et al. 3D Additive Manufactured Composite Scaffolds with Antibiotic-Loaded Lamellar Fillers for Bone Infection Prevention and Tissue Regeneration. Bioact. Mater. 2021, 6, 1073–1082. [Google Scholar] [CrossRef]
  41. Bastianini, M.; Scatto, M.; Sisani, M.; Scopece, P.; Patelli, A.; Petracci, A. Innovative Composites Based on Organic Modified Zirconium Phosphate and PEOT/PBT Copolymer. J. Compos. Sci. 2018, 2, 31. [Google Scholar] [CrossRef] [Green Version]
  42. F04 Committee. ASTM F2027-08, Standard Guide for Characterization and Testing of Raw or Starting Biomaterials for Tissue-Engineered Medical Products; ASTM International: West Conshohocken, PA, USA, 2008. [Google Scholar]
  43. D20 Committee. ASTM D638-08, Standard Test Method for Tensile Properties of Plastics; ASTM International: West Conshohocken, PA, USA, 2008. [Google Scholar]
  44. D20 Committee. ASTM D695-10, Standard Test Method for Compressive Properties of Rigid Plastics; ASTM International: West Conshohocken, PA, USA, 2010. [Google Scholar]
  45. Technical Committee: ISO/TC 61/SC 2 Mechanical Behavior. ISO 527-1:2019 Plastics—Determination of Tensile Properties—Part 1: General Principles; IOS: Geneva, Switzerland, 2019; p. 26. [Google Scholar]
  46. Technical Committee: ISO/TC 61/SC 2 Mechanical Behavior. ISO 527-2:2012 Plastics—Determination of Tensile Properties—Part 2: Test Conditions for Moulding and Extrusion Plastics; IOS: Geneva, Switzerland, 2012; p. 11. [Google Scholar]
  47. Technical Committee: ISO/TC 61/SC 13 Composites and Reinforcement Fibres. ISO 527-4:1997 Plastics—Determination of Tensile Properties—Part 4: Test Conditions for Isotropic and Orthotropic Fibre-Reinforced Plastic Composites; IOS: Geneva, Switzerland, 1997; p. 11. [Google Scholar]
  48. F04 Committee. ASTM F2150-02e1, Standard Guide for Characterization and Testing of Biomaterial Scaffolds Used in Tissue-Engineered Medical Products; ASTM International: West Conshohocken, PA, USA, 2002. [Google Scholar]
  49. Technical Committee: ISO/TC 61/SC 2 Mechanical Behavior. ISO 604:2002 Plastics—Determination of Compressive Properties; IOS: Geneva, Switzerland, 2002; p. 18. [Google Scholar]
  50. Calderón-Uríszar-Aldaca, I.; Biezma, M.V.; Matanza, A.; Briz, E.; Bastidas, D.M. Second-Order Fatigue of Intrinsic Mean Stress under Random Loadings. Int. J. Fatigue 2020, 130, 105257. [Google Scholar] [CrossRef]
  51. Calderon-Uriszar-Aldaca, I.; Briz, E.; Biezma, M.V.; Puente, I. A Plain Linear Rule for Fatigue Analysis under Natural Loading Considering the Coupled Fatigue and Corrosion Effect. Int. J. Fatigue 2019, 122, 141–151. [Google Scholar] [CrossRef]
  52. Calderon-Uriszar-Aldaca, I.; Biezma, M.V. A Plain Linear Rule for Fatigue Analysis under Natural Loading Considering the Sequence Effect. Int. J. Fatigue 2017, 103, 386–394. [Google Scholar] [CrossRef] [Green Version]
Figure 2. Finite element modelling of 3D-printed bone scaffolds with variable properties depending on bone topology: (A) Determination of equivalent finite element and bone discretization; (B) Flowchart of required steps for realistic finite element modelling of bone scaffolds within a bone.
Figure 2. Finite element modelling of 3D-printed bone scaffolds with variable properties depending on bone topology: (A) Determination of equivalent finite element and bone discretization; (B) Flowchart of required steps for realistic finite element modelling of bone scaffolds within a bone.
Mathematics 09 01746 g002
Figure 3. (a) Compression samples of PEOT/PBT 45% nano HA; (b) Tension samples made of different nano HA fillers concentration in PEOT/PBT.
Figure 3. (a) Compression samples of PEOT/PBT 45% nano HA; (b) Tension samples made of different nano HA fillers concentration in PEOT/PBT.
Mathematics 09 01746 g003
Figure 4. Compression stress–strain curve of PEOT/PBT with 45% nanoHA samples and mean behaviour considered for FEM.
Figure 4. Compression stress–strain curve of PEOT/PBT with 45% nanoHA samples and mean behaviour considered for FEM.
Mathematics 09 01746 g004
Figure 5. Standardized compression tests over 3D printed scaffold made of PEOT/PBT 45% nano HA.
Figure 5. Standardized compression tests over 3D printed scaffold made of PEOT/PBT 45% nano HA.
Mathematics 09 01746 g005
Figure 6. Curves of Stress (MPa)—Strain (%) relationship of three scaffold types (a) G25_nanoHA 1.25; (b) G25_nanoHA_1.00; (c) G25_nanoHA_0.75.
Figure 6. Curves of Stress (MPa)—Strain (%) relationship of three scaffold types (a) G25_nanoHA 1.25; (b) G25_nanoHA_1.00; (c) G25_nanoHA_0.75.
Mathematics 09 01746 g006
Figure 7. Geometry of scaffold numerical models: (a) Global view showing double symmetry of G25_nanoHA 1.25; (b) Cross-section of the resulting scaffold G25_nanoHA 1.25; (c) Global view showing double symmetry G25_nanoHA_1.00; (d) Cross-section of the resulting scaffold G25_nanoHA_1.00; (e) Global view showing double symmetry G25_nanoHA_1.00; (f) Cross-section of the resulting scaffold G25_nanoHA_0.75.
Figure 7. Geometry of scaffold numerical models: (a) Global view showing double symmetry of G25_nanoHA 1.25; (b) Cross-section of the resulting scaffold G25_nanoHA 1.25; (c) Global view showing double symmetry G25_nanoHA_1.00; (d) Cross-section of the resulting scaffold G25_nanoHA_1.00; (e) Global view showing double symmetry G25_nanoHA_1.00; (f) Cross-section of the resulting scaffold G25_nanoHA_0.75.
Mathematics 09 01746 g007
Figure 8. Detail of the meshing for G25_nanoHA_1.25 finite element model.
Figure 8. Detail of the meshing for G25_nanoHA_1.25 finite element model.
Mathematics 09 01746 g008
Figure 9. Boundary conditions considered for FEM, only G25_nanoHA_0.75 case is shown: (a) Detail of bonding between perpendicular fibres; (b) Detail of the lower; (c) Detail of the upper plane.
Figure 9. Boundary conditions considered for FEM, only G25_nanoHA_0.75 case is shown: (a) Detail of bonding between perpendicular fibres; (b) Detail of the lower; (c) Detail of the upper plane.
Mathematics 09 01746 g009
Figure 10. Results of FEM simulation of a monotonic compression test on a G25_nanoHA_1.25 Scaffold: (a) Stress–strain curve; (b) Nodal stresses (Von Mises equivalent stress, in Pascals).
Figure 10. Results of FEM simulation of a monotonic compression test on a G25_nanoHA_1.25 Scaffold: (a) Stress–strain curve; (b) Nodal stresses (Von Mises equivalent stress, in Pascals).
Mathematics 09 01746 g010
Figure 11. Results of FEM simulation of a monotonic compression test on a G25_nanoHA_1.00 Scaffold: (a) Stress–strain curve; (b) Nodal stresses (Von Mises equivalent stress, in 10−1 Pascals).
Figure 11. Results of FEM simulation of a monotonic compression test on a G25_nanoHA_1.00 Scaffold: (a) Stress–strain curve; (b) Nodal stresses (Von Mises equivalent stress, in 10−1 Pascals).
Mathematics 09 01746 g011
Figure 12. Results of FEM simulation of a monotonic compression test on a G25_nanoHA_0.75 Scaffold: (a) Stress–strain curve; (b) Nodal stresses (Von Mises equivalent stress, in 10−1 Pascals).
Figure 12. Results of FEM simulation of a monotonic compression test on a G25_nanoHA_0.75 Scaffold: (a) Stress–strain curve; (b) Nodal stresses (Von Mises equivalent stress, in 10−1 Pascals).
Mathematics 09 01746 g012
Figure 13. Influence of the Poisson’s ratio on the stress–strain relationships from FEM.
Figure 13. Influence of the Poisson’s ratio on the stress–strain relationships from FEM.
Mathematics 09 01746 g013
Figure 14. Configurations considered to analyse the influence of the position of the fibres on the compressive stiffness of the scaffold (G25_nanoHA 1.25).
Figure 14. Configurations considered to analyse the influence of the position of the fibres on the compressive stiffness of the scaffold (G25_nanoHA 1.25).
Mathematics 09 01746 g014
Figure 15. Geometry of the numerical models for configurations 1 and 2: (a) Global view of configuration 1; (b) Cross-section of configuration 1; (c) Global view configuration 2; (d) Cross-section of configuration 2.
Figure 15. Geometry of the numerical models for configurations 1 and 2: (a) Global view of configuration 1; (b) Cross-section of configuration 1; (c) Global view configuration 2; (d) Cross-section of configuration 2.
Mathematics 09 01746 g015
Figure 16. Influence of the relative position of the fibres within the scaffold G25_nanoHA_1.25 on the stress–strain curve.
Figure 16. Influence of the relative position of the fibres within the scaffold G25_nanoHA_1.25 on the stress–strain curve.
Mathematics 09 01746 g016
Figure 17. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) Configuration 1, in Pascals; (b) Configuration 2, in 10−1 Pascals.
Figure 17. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) Configuration 1, in Pascals; (b) Configuration 2, in 10−1 Pascals.
Mathematics 09 01746 g017
Figure 18. Influence of the fibre diameter on the stress–strain relationships from FEM.
Figure 18. Influence of the fibre diameter on the stress–strain relationships from FEM.
Mathematics 09 01746 g018
Figure 19. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) φfibre = 0.235 mm, in 10−1 Pascals; (b) φfibre = 0.220 mm, in 10−1 Pascals.
Figure 19. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) φfibre = 0.235 mm, in 10−1 Pascals; (b) φfibre = 0.220 mm, in 10−1 Pascals.
Mathematics 09 01746 g019
Figure 20. Influence of the strand distance on the stress–strain relationships from FEM.
Figure 20. Influence of the strand distance on the stress–strain relationships from FEM.
Mathematics 09 01746 g020
Figure 21. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) d = 1.20 mm, in 10−1 Pascals; (b) d = 1.30 mm, in 10−1 Pascals.
Figure 21. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) d = 1.20 mm, in 10−1 Pascals; (b) d = 1.30 mm, in 10−1 Pascals.
Mathematics 09 01746 g021
Figure 22. Influence of the layer height on the stress–strain relationships from FEM.
Figure 22. Influence of the layer height on the stress–strain relationships from FEM.
Mathematics 09 01746 g022
Figure 23. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) h = 0.215 mm, in 10−1 Pascals; (b) h = 0.230 mm, in 10−1 Pascals.
Figure 23. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) h = 0.215 mm, in 10−1 Pascals; (b) h = 0.230 mm, in 10−1 Pascals.
Mathematics 09 01746 g023
Figure 24. Influence of the sample diameter on the stress–strain relationships from FEM.
Figure 24. Influence of the sample diameter on the stress–strain relationships from FEM.
Mathematics 09 01746 g024
Figure 25. Detail of printed scaffolds.
Figure 25. Detail of printed scaffolds.
Mathematics 09 01746 g025
Figure 26. Influence of the deflection of the printed fibres on the stress–strain relationships from FEM.
Figure 26. Influence of the deflection of the printed fibres on the stress–strain relationships from FEM.
Mathematics 09 01746 g026
Figure 27. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) f = 0.06 mm, in 10−1 Pascals; (b) f = 0.12 mm, in 10−1 Pascals.
Figure 27. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): (a) f = 0.06 mm, in 10−1 Pascals; (b) f = 0.12 mm, in 10−1 Pascals.
Mathematics 09 01746 g027
Figure 28. Location of fibres in situation 1.
Figure 28. Location of fibres in situation 1.
Mathematics 09 01746 g028
Figure 29. Geometry of the numerical model for situation 1 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section.
Figure 29. Geometry of the numerical model for situation 1 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section.
Mathematics 09 01746 g029
Figure 30. Location of fibres in situation 2.
Figure 30. Location of fibres in situation 2.
Mathematics 09 01746 g030
Figure 31. Geometry of the numerical model for situation 2 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section.
Figure 31. Geometry of the numerical model for situation 2 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section.
Mathematics 09 01746 g031
Figure 32. Influence of the straightness of the columns on the elastic modulus from FEM of G25_nanoHA 1.25: (a) Influence on the stress–strain relationships of configuration 1; (b) Influence on the stress–strain relationships of configuration 2.
Figure 32. Influence of the straightness of the columns on the elastic modulus from FEM of G25_nanoHA 1.25: (a) Influence on the stress–strain relationships of configuration 1; (b) Influence on the stress–strain relationships of configuration 2.
Mathematics 09 01746 g032aMathematics 09 01746 g032b
Figure 33. Influence of the existence of broken fibres on the stress–strain curve (G25_nanoHA 1.25).
Figure 33. Influence of the existence of broken fibres on the stress–strain curve (G25_nanoHA 1.25).
Mathematics 09 01746 g033
Figure 34. Gaussian density function of each uncertainty reduction factor.
Figure 34. Gaussian density function of each uncertainty reduction factor.
Mathematics 09 01746 g034
Figure 35. Gaussian density function of global uncertainty factor derived from variable product.
Figure 35. Gaussian density function of global uncertainty factor derived from variable product.
Mathematics 09 01746 g035
Table 1. Geometrical properties of the tested scaffolds to be represented in the numerical models.
Table 1. Geometrical properties of the tested scaffolds to be represented in the numerical models.
Scaffold TypeFibre Diameter
(mm)
Strand Distance
(mm)
Layer Height
(mm)
Scaffold Diameter
(mm)]
Scaffold Height
(mm)
G25_nanoHA 1.250.251.250.204.004.05
G25_nanoHA_1.000.251.000.204.004.05
G25_nanoHA_0.750.250.750.204.004.05
Table 2. Results from monotonic compression tests of each sample and type, with corresponding mean value (Mean) and standard deviation (SD) of each type.
Table 2. Results from monotonic compression tests of each sample and type, with corresponding mean value (Mean) and standard deviation (SD) of each type.
SampleYoung Modulus
(MPa)
Yield Strength
(MPa)
Elongation at Yield
(%)
G25_nanoHA 1.25_12.330.259.80
G25_nanoHA 1.25_21.860.2512.20
G25_nanoHA 1.25_32.560.2911.60
G25_nanoHA_1.25 (Mean)2.250.2611.20
G25_nanoHA_1.25 (SD)0.360.021.25
G25_nanoHA_1.00_15.540.4210.90
G25_nanoHA_1.00_23.980.5216.00
G25_nanoHA_1.00_35.880.5013.50
G25_nanoHA_1.00 (Mean)5.130.4813.47
G25_nanoHA_1.00 (SD)1.010.052.55
G25_nanoHA_0.75_16.930.7815.90
G25_nanoHA_0.75_26.170.8016.90
G25_nanoHA_0.75_35.900.7916.40
G25_nanoHA_0.75 (Mean)6.330.7916.40
G25_nanoHA_0.75 (SD)0.530.010.50
Table 3. Comparison between the Young’s modulus obtained using FEM and direct scaffold testing, with corresponding Mean and Standard Deviation values.
Table 3. Comparison between the Young’s modulus obtained using FEM and direct scaffold testing, with corresponding Mean and Standard Deviation values.
ScaffoldG25 nanoHA 1.25G25 nanoHA 1.00G25 nanoHA 0.75
ModelMeanSDModelMeanSDModelMeanSD
Young’s Modulus(MPa)6.602.250.367.905.131.0115.496.330.53
Table 4. Comparison between base G25_nanoHA_1.25 characteristics and developed variations for sensitivity analysis.
Table 4. Comparison between base G25_nanoHA_1.25 characteristics and developed variations for sensitivity analysis.
SourceBaseVariation 1Variation 2
Poisson’s Ratio0.40.3750.35
Position of fibresAt plane of sym.Out plane of sym.-
Fibre diameter (mm)0.250.2350.22
Strand distance (mm)1.251.21.3
Layer height (mm)0.20.2150.23
Sample diameter (mm)43.94.1
Deflection of fibres (mm)00.060.12
Straightness of columnsStraightAlternatedCumulated
Existence of broken fibresNo broken fibresbreakage at crossingfar from crossing
Table 5. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the Poisson’s ratio.
Table 5. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the Poisson’s ratio.
Poisson’s Ratio (-)Elastic Modulus (MPa)
0.4006.60
0.3756.57
0.3506.54
Table 6. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different configurations in terms of the relative position of the fibres within the scaffold.
Table 6. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different configurations in terms of the relative position of the fibres within the scaffold.
Poisson’s Ratio (-)Elastic Modulus (MPa)
0.4006.60
0.3756.26
Table 7. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the fibre diameter.
Table 7. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the fibre diameter.
Fibre Diameter (mm)Elastic Modulus (MPa)
0.2506.60
0.2355.08
0.2203.39
Table 8. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the layer height.
Table 8. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the layer height.
Strand Distance (mm)Elastic Modulus (MPa)
1.206.61
1.256.60
1.306.58
Table 9. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the layer height.
Table 9. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the layer height.
Layer Height (mm)Elastic Modulus (MPa)
0.2006.60
0.2155.53
0.2304.07
Table 10. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the sample diameter. The theoretical value for the sample diameter has been used for calculating the stress.
Table 10. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the sample diameter. The theoretical value for the sample diameter has been used for calculating the stress.
Sample Diameter (mm)Elastic Modulus (MPa)
3.906.61
4.006.60
4.106.62
Table 11. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the deflection of the printed fibres.
Table 11. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the deflection of the printed fibres.
Deflection of the Printed Fibres (mm)Elastic Modulus (MPa)
0.006.60
0.066.45
0.126.41
Table 12. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the eccentricity or deviation in situation 1.
Table 12. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the eccentricity or deviation in situation 1.
Deviation (mm)Elastic Modulus (MPa)
0.0006.70
0.0255.76
0.0503.99
Table 13. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the eccentricity or deviation in situation 2.
Table 13. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the eccentricity or deviation in situation 2.
Deviation (mm)Elastic Modulus (MPa)
0.0006.70
0.0056.69
0.0096.67
Table 14. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different situations concerning the existence of broken fibres.
Table 14. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different situations concerning the existence of broken fibres.
Situation (#)Elastic Modulus (MPa)
Situation 16.70
Situation 26.57
Situation 36.43
Situation 46.70
Situation 56.70
Table 15. Uncertainty factors Fi and corresponding mean μi and standard deviations σi, whose minimum and maximum values correspond to percentiles 5 and 95, respectively.
Table 15. Uncertainty factors Fi and corresponding mean μi and standard deviations σi, whose minimum and maximum values correspond to percentiles 5 and 95, respectively.
Variable
(xi)
Source
(-)
Fi,min
(#)
Fi,max
(#)
μi
(#)
σi
(#)
x1Poisson’s Ratio0.99091.00000.99550.0028
x2Position of the fibres0.94851.00000.97420.0157
x3Fibre Diameter0.51361.00000.75680.1478
x4Strand Distance0.99551.00000.99770.0014
x5Layer Height0.61671.00000.80830.1165
x6Sample Diameter0.95451.05301.00380.0299
x7Deflection of the fibres0.97121.00000.98560.0088
x8Straightness of the columns0.59551.00000.79780.1230
x9Existence of broken fibres0.95971.00000.97990.0122
Table 16. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone without correction factor (γMP = 1).
Table 16. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone without correction factor (γMP = 1).
EA,FEMEm,dεmδEFEMσFEM
(MPa)(MPa)(-)(mm)(MPa)(MPa)
G25 nanoHA 1.256.66.60.0591.182566.233.4
G25 nanoHA 1.007.97.90.0490.987730.736.1
G25 nanoHA 0.7515.4915.50.0250.503563.114.2
Table 17. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with mean correction factor (γMP = 2.2).
Table 17. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with mean correction factor (γMP = 2.2).
EA,FEMEm,dεmδEFEMσFEM
(MPa)(MPa)(-)(mm)(MPa)(MPa)
G25 nanoHA 1.256.63.00.1302.600566.273.6
G25 nanoHA 1.007.93.60.1092.172730.779.3
G25 nanoHA 0.7515.497.00.0551.108563.131.2
Table 18. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with characteristic correction factor (γMP = 4.25).
Table 18. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with characteristic correction factor (γMP = 4.25).
EA,FEMEk,dεkδEFEMσFEM
(MPa)(MPa)(-)(mm)(MPa)(MPa)
G25 nanoHA 1.256.61.60.2515.022566.2142.2
G25 nanoHA 1.007.91.90.2104.195730.7153.3
G25 nanoHA 0.7515.493.60.1072.140563.160.2
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Calderon-Uriszar-Aldaca, I.; Perez, S.; Sinha, R.; Camara-Torres, M.; Villanueva, S.; Mota, C.; Patelli, A.; Matanza, A.; Moroni, L.; Sanchez, A. Size Effects in Finite Element Modelling of 3D Printed Bone Scaffolds Using Hydroxyapatite PEOT/PBT Composites. Mathematics 2021, 9, 1746. https://0-doi-org.brum.beds.ac.uk/10.3390/math9151746

AMA Style

Calderon-Uriszar-Aldaca I, Perez S, Sinha R, Camara-Torres M, Villanueva S, Mota C, Patelli A, Matanza A, Moroni L, Sanchez A. Size Effects in Finite Element Modelling of 3D Printed Bone Scaffolds Using Hydroxyapatite PEOT/PBT Composites. Mathematics. 2021; 9(15):1746. https://0-doi-org.brum.beds.ac.uk/10.3390/math9151746

Chicago/Turabian Style

Calderon-Uriszar-Aldaca, Iñigo, Sergio Perez, Ravi Sinha, Maria Camara-Torres, Sara Villanueva, Carlos Mota, Alessandro Patelli, Amaia Matanza, Lorenzo Moroni, and Alberto Sanchez. 2021. "Size Effects in Finite Element Modelling of 3D Printed Bone Scaffolds Using Hydroxyapatite PEOT/PBT Composites" Mathematics 9, no. 15: 1746. https://0-doi-org.brum.beds.ac.uk/10.3390/math9151746

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop