Size effects in finite element modelling of 3D printed bone scaffolds using hydroxyapatite PEOT/PBT composites

mathematics Article Size Effects in Finite Element Modelling of 3D Printed Bone Scaffolds Using Hydroxyapatite PEOT/PBT Composites Iñigo Calderon-Uriszar-Aldaca 1,2,*, Sergio Perez 1, Ravi Sinha 3, Maria Camara-Torres 3, Sara Villanueva 1, Carlos Mota 3, Alessandro Patelli 4, Amaia Matanza 5, Lorenzo Moroni 3and Alberto Sanchez 1,*   Citation: 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://doi.org/ 10.3390/math9151746 Academic Editor: Mauro Malvè Received: 2 June 2021 Accepted: 7 July 2021 Published: 24 July 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). 1TECNALIA, Basque Research and Technology Alliance (BRTA), Mikeletegi Pasealekua 2, 20009 Donostia-San Sebastián, Spain; sergio.per[email protected] (S.P.); [email protected] (S.V.) 2Department 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; [email protected] (R.S.); [email protected] (M.C.-T.); [email protected] (C.M.); [email protected] (L.M.) 4Department of Physics and Astronomy, Padova University, Via Marzolo 8, 35131 Padova, Italy; [email protected] 5Materials Physics Center (MPC), Centro de Fisica de Materiales (CSIC, UPV/EHU), Paseo Manuel de Lardizabal 5, 20018 Donostia-San Sebastián, Spain; [email protected] *Correspondence: [email protected] (I.C.-U.-A.); [email protected] (A.S.) Abstract: Additive manufacturing (AM) of scaffolds enables the fabrication of customized patientspecific 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 cm 3 ). 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 cm 3 , 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. Keywords: finite element modelling; bone tissue engineering; 3D scaffold; additive manufacturing 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]. Mathematics 2021,9, 1746. https://doi.org/10.3390/math9151746 https://www.mdpi.com/journal/mathematics Mathematics 2021,9, 1746 2 of 37 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 – 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 – 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 cm 3 ). 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. Mathematics 2021, 9, 1746 2 of 37 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–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–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. (a) (b) 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 2021,9, 1746 3 of 37 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 microscale 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 – 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). Mathematics 2021,9, 1746 4 of 37 Mathematics 2021, 9, 1746 4 of 37 (A) (B) 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. 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–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. 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. 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 – 47 ] (Figure 4). A Poisson’s ratio of 0.4 was taken as usual practice for polymers, although Mathematics 2021,9, 1746 5 of 37 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. Mathematics 2021, 9, 1746 5 of 37 (a) (b) 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 4. Compression stress–strain curve of PEOT/PBT with 45% nanoHA samples and mean behaviour considered for FEM. 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). Figure 5. Standardized compression tests over 3D printed scaffold made of PEOT/PBT 45% nano HA. 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 2021, 9, 1746 5 of 37 (a) (b) 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 4. Compression stress–strain curve of PEOT/PBT with 45% nanoHA samples and mean behaviour considered for FEM. 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). Figure 5. Standardized compression tests over 3D printed scaffold made of PEOT/PBT 45% nano HA. Figure 4. Compression stress–strain curve of PEOT/PBT with 45% nanoHA samples and mean behaviour considered for FEM. 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). Mathematics 2021,9, 1746 6 of 37 Mathematics 2021, 9, 1746 5 of 37 (a) (b) 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 4. Compression stress–strain curve of PEOT/PBT with 45% nanoHA samples and mean behaviour considered for FEM. 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). 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. 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. Table 1. Geometrical properties of the tested scaffolds to be represented in the numerical models. Scaffold Type Fibre Diameter (mm) Strand Distance (mm) Layer Height (mm) Scaffold Diameter (mm)] Scaffold Height (mm) G25_nanoHA 1.25 0.25 1.25 0.20 4.00 4.05 G25_nanoHA_1.00 0.25 1.00 0.20 4.00 4.05 G25_nanoHA_0.75 0.25 0.75 0.20 4.00 4.05 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. Mathematics 2021,9, 1746 7 of 37 Mathematics 2021, 9, 1746 6 of 37 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. Table 1. Geometrical properties of the tested scaffolds to be represented in the numerical models. Scaffold Type Fibre Diameter (mm) Strand Distance (mm) Layer Height (mm) Scaffold Diameter (mm)] Scaffold Height (mm) G25_nanoHA 1.25 0.25 1.25 0.20 4.00 4.05 G25_nanoHA_1.00 0.25 1.00 0.20 4.00 4.05 G25_nanoHA_0.75 0.25 0.75 0.20 4.00 4.05 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. (a) (b) (c) 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. 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. Table 2. Results from monotonic compression tests of each sample and type, with corresponding mean value (Mean) and standard deviation (SD) of each type. Sample Young Modulus (MPa) Yield Strength (MPa) Elongation at Yield (%) G25_nanoHA 1.25_1 2.33 0.25 9.80 G25_nanoHA 1.25_2 1.86 0.25 12.20 G25_nanoHA 1.25_3 2.56 0.29 11.60 G25_nanoHA_1.25 (Mean) 2.25 0.26 11.20 G25_nanoHA_1.25 (SD) 0.36 0.02 1.25 G25_nanoHA_1.00_1 5.54 0.42 10.90 G25_nanoHA_1.00_2 3.98 0.52 16.00 G25_nanoHA_1.00_3 5.88 0.50 13.50 G25_nanoHA_1.00 (Mean) 5.13 0.48 13.47 G25_nanoHA_1.00 (SD) 1.01 0.05 2.55 G25_nanoHA_0.75_1 6.93 0.78 15.90 G25_nanoHA_0.75_2 6.17 0.80 16.90 G25_nanoHA_0.75_3 5.90 0.79 16.40 G25_nanoHA_0.75 (Mean) 6.33 0.79 16.40 G25_nanoHA_0.75 (SD) 0.53 0.01 0.50 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 Mathematics 2021,9, 1746 8 of 37 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 2in 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. Mathematics 2021, 9, 1746 8 of 37 (a) (b) (c) (d) (e) (f) 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. 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. 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. 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 Mathematics 2021,9, 1746 9 of 37 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. Mathematics 2021, 9, 1746 9 of 37 Figure 8. Detail of the meshing for G25_nanoHA_1.25 finite element model. 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). (a) (b) (c) 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 8. Detail of the meshing for G25_nanoHA_1.25 finite element model. 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. Mathematics 2021,9, 1746 16 of 37 Mathematics 2021, 9, 1746 15 of 37 (a) (b) (c) (d) 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. 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). Figure 16. Influence of the relative position of the fibres within the scaffold G25_nanoHA_1.25 on the stress–strain curve. 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). Figure 16. Influence of the relative position of the fibres within the scaffold G25_nanoHA_1.25 on the stress–strain curve. 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). 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.400 6.60 0.375 6.26 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. Mathematics 2021, 9, 1746 16 of 37 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.400 6.60 0.375 6.26 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. (a) (b) 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. 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). Figure 17. Nodal stresses of G25_nanoHA_1.25 (Von Mises equivalent stress): ( a ) Configuration 1, in Pascals; ( b ) Configuration 2, in 10−1Pascals. Mathematics 2021,9, 1746 17 of 37 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). Mathematics 2021, 9, 1746 16 of 37 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.400 6.60 0.375 6.26 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. (a) (b) 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. 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). Figure 18. Influence of the fibre diameter on the stress–strain relationships from FEM. 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. 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.250 6.60 0.235 5.08 0.220 3.39 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 7and 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). Mathematics 2021,9, 1746 18 of 37 Mathematics 2021, 9, 1746 17 of 37 Figure 18. Influence of the fibre diameter on the stress–strain relationships from FEM. 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. 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.250 6.60 0.235 5.08 0.220 3.39 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). (a) (b) 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. 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. 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−1Pascals. 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. Mathematics 2021, 9, 1746 18 of 37 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.20 6.61 1.25 6.60 1.30 6.58 Figure 20. Influence of the strand distance on the stress–strain relationships from FEM. 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). (a) (b) 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 20. Influence of the strand distance on the stress–strain relationships from FEM. 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. 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.20 6.61 1.25 6.60 1.30 6.58 Mathematics 2021,9, 1746 19 of 37 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). Mathematics 2021, 9, 1746 18 of 37 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.20 6.61 1.25 6.60 1.30 6.58 Figure 20. Influence of the strand distance on the stress–strain relationships from FEM. 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). (a) (b) 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−1Pascals. 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). Mathematics 2021, 9, 1746 19 of 37 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). Figure 22. Influence of the layer height on the stress–strain relationships from FEM. 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. 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.200 6.60 0.215 5.53 0.230 4.07 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). Figure 22. Influence of the layer height on the stress–strain relationships from FEM. Mathematics 2021,9, 1746 20 of 37 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. 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.200 6.60 0.215 5.53 0.230 4.07 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 9and 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). Mathematics 2021, 9, 1746 20 of 37 (a) (b) 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. 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). Figure 24. Influence of the sample diameter on the stress–strain relationships from FEM. 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. 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−1Pascals. 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). Mathematics 2021,9, 1746 21 of 37 Mathematics 2021, 9, 1746 20 of 37 (a) (b) 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. 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). Figure 24. Influence of the sample diameter on the stress–strain relationships from FEM. 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. Figure 24. Influence of the sample diameter on the stress–strain relationships from FEM. 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. 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.90 6.61 4.00 6.60 4.10 6.62 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. Mathematics 2021,9, 1746 22 of 37 Mathematics 2021, 9, 1746 21 of 37 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.90 6.61 4.00 6.60 4.10 6.62 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. Figure 25. Detail of printed scaffolds. 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). Figure 25. Detail of printed scaffolds. 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). Mathematics 2021, 9, 1746 22 of 37 Figure 26. Influence of the deflection of the printed fibres on the stress–strain relationships from FEM. 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. 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.00 6.60 0.06 6.45 0.12 6.41 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. Figure 26. Influence of the deflection of the printed fibres on the stress–strain relationships from FEM. 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. 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.00 6.60 0.06 6.45 0.12 6.41 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 Mathematics 2021,9, 1746 23 of 37 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. Mathematics 2021, 9, 1746 23 of 37 (a) (b) 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. 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. 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−1Pascals. 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. Mathematics 2021,9, 1746 24 of 37 Mathematics 2021, 9, 1746 24 of 37 Figure 28. Location of fibres in situation 1. 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. (a) (b) Figure 29. Geometry of the numerical model for situation 1 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section. 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. Figure 30. Location of fibres in situation 2. Figure 28. Location of fibres in situation 1. 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. Mathematics 2021, 9, 1746 24 of 37 Figure 28. Location of fibres in situation 1. 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. (a) (b) Figure 29. Geometry of the numerical model for situation 1 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section. 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. Figure 30. Location of fibres in situation 2. Figure 29. Geometry of the numerical model for situation 1 for G25_nanoHA 1.25: ( a ) Global view; ( b ) Vertical plane section. 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. Mathematics 2021, 9, 1746 24 of 37 Figure 28. Location of fibres in situation 1. 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. (a) (b) Figure 29. Geometry of the numerical model for situation 1 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section. 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. Figure 30. Location of fibres in situation 2. Figure 30. Location of fibres in situation 2. Mathematics 2021,9, 1746 25 of 37 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. Mathematics 2021, 9, 1746 25 of 37 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. (a) (b) Figure 31. Geometry of the numerical model for situation 2 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section. 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 Figures 29 and 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). (a) Figure 31. Geometry of the numerical model for situation 2 for G25_nanoHA 1.25: ( a ) Global view; ( b ) Vertical plane section. 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 Figures 29 and 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). Mathematics 2021, 9, 1746 25 of 37 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. (a) (b) Figure 31. Geometry of the numerical model for situation 2 for G25_nanoHA 1.25: (a) Global view; (b) Vertical plane section. 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 Figures 29 and 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). (a) Figure 32. Cont. Mathematics 2021,9, 1746 32 of 37 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[ˆ x1·ˆ x2··· ˆ x9]≈0.4578 (19) σ=qVAR(x1·x2·x3··· x9)≈0.1346 (20) 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. Mathematics 2021, 9, 1746 32 of 37 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. Figure 35. Gaussian density function of global uncertainty factor derived from variable product. 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: Ed=EM γMP (21) Ed=0.4578∙EM≈EM 2.2 (22) Ed=0.2364∙EM≈EM 4.25 (23) 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. 0 0.5 1 1.5 2 2.5 3 3.5 0 0.06 0.12 0.18 0.24 0.3 0.36 0.42 0.48 0.54 0.6 0.66 0.72 0.78 0.84 0.9 0.96 1.02 1.08 1.14 1.2 f(x1∙x2∙…∙x9) X Figure 35. Gaussian density function of global uncertainty factor derived from variable product. 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 E d can be derived from the elastic modulus obtained from models with idealized or theoretical conditions, E M , 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: Ed=EM γMP (21) Ed=0.4578·EM≈EM 2.2 (22) Ed=0.2364·EM≈EM 4.25 (23) Mathematics 2021,9, 1746 33 of 37 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 mm 2 . 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 Ex−µyx·σy Ey−µzx·σz Ez =σx Ex(24) Now, using the apparent young modulus obtained from the idealized models E A,FEM already presented in Table 3without any correction factor, the mean young modulus, E m,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, E FEM , 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. Table 16. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone without correction factor (γMP = 1). EA,FEM Em,d εmδEFEM σFEM (MPa) (MPa) (-) (mm) (MPa) (MPa) G25 nanoHA 1.25 6.6 6.6 0.059 1.182 566.2 33.4 G25 nanoHA 1.00 7.9 7.9 0.049 0.987 730.7 36.1 G25 nanoHA 0.75 15.49 15.5 0.025 0.503 563.1 14.2 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 E m,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. Mathematics 2021,9, 1746 34 of 37 Table 17. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with mean correction factor (γMP = 2.2). EA,FEM Em,d εmδEFEM σFEM (MPa) (MPa) (-) (mm) (MPa) (MPa) G25 nanoHA 1.25 6.6 3.0 0.130 2.600 566.2 73.6 G25 nanoHA 1.00 7.9 3.6 0.109 2.172 730.7 79.3 G25 nanoHA 0.75 15.49 7.0 0.055 1.108 563.1 31.2 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 E k,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). Table 18. Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with characteristic correction factor (γMP = 4.25). EA,FEM Ek,d εkδEFEM σFEM (MPa) (MPa) (-) (mm) (MPa) (MPa) G25 nanoHA 1.25 6.6 1.6 0.251 5.022 566.2 142.2 G25 nanoHA 1.00 7.9 1.9 0.210 4.195 730.7 153.3 G25 nanoHA 0.75 15.49 3.6 0.107 2.140 563.1 60.2 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. Mathematics 2021,9, 1746 35 of 37 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 cm 3 ) 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. Mathematics 2021,9, 1746 36 of 37 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. [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. [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. [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. [CrossRef] 5. Karageorgiou, V.; Kaplan, D. Porosity of 3D Biomaterial Scaffolds and Osteogenesis. Biomaterials 2005 ,26, 5474–5491. [CrossRef] 6. Hutmacher, D.W. Scaffolds in Tissue Engineering Bone and Cartilage. Biomaterials 2000,21, 2529–2543. [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. [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. [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. [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. [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. [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. [CrossRef] 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. [CrossRef] 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. [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. [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. [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. [CrossRef] [PubMed] 18. Tjong, S.C. Structural and Mechanical Properties of Polymer Nanocomposites. Mater. Sci. Eng. R Rep. 2006 ,53, 73–197. [CrossRef] 19. Paul, D.R.; Robeson, L.M. Polymer Nanotechnology: Nanocomposites. Polymer 2008,49, 3187–3204. [CrossRef] 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. [CrossRef] 21. Pavlidou, S.; Papaspyrides, C.D. A Review on Polymer–Layered Silicate Nanocomposites. Prog. Polym. Sci. 2008 ,33, 1119–1198. [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. [CrossRef] 23. Potts, J.R.; Dreyer, D.R.; Bielawski, C.W.; Ruoff, R.S. Graphene-Based Polymer Nanocomposites. Polymer 2011 ,52, 5–25. [CrossRef] 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. [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. [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. [CrossRef] Mathematics 2021,9, 1746 37 of 37 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. [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. [CrossRef] 29. Munarin, F.; Petrini, P.; Gentilini, R.; Pillai, R.S.; Dirè, S.; Tanzi, M.C.; Sglavo, V.M. Microand Nano-Hydroxyapatite as Active Reinforcement for Soft Biocomposites. Int. J. Biol. Macromol. 2015,72, 199–209. [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. [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. [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. [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. [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. [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. [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. [CrossRef] 37. Almeida, H.A. Design of Tissue Engineering Scaffolds Based on Hyperbolic Surfaces: Structural Numerical Evaluation. Med. Eng. 2014,8, 1033–1040. [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. [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. [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. [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. [CrossRef] 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. 43. D20 Committee. ASTM D638-08, Standard Test Method for Tensile Properties of Plastics; ASTM International: West Conshohocken, PA, USA, 2008. 44. D20 Committee. ASTM D695-10, Standard Test Method for Compressive Properties of Rigid Plastics; ASTM International: West Conshohocken, PA, USA, 2010. 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. 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. 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. 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. 49. Technical Committee: ISO/TC 61/SC 2 Mechanical Behavior. ISO 604:2002 Plastics—Determination of Compressive Properties; IOS: Geneva, Switzerland, 2002; p. 18. 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. [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. [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. [CrossRef]