*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 <sup>3</sup> and Alberto Sanchez 1,\***


**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].

Academic Editor: Mauro Malvè

Received: 2 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.


**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

**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/).

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.

**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].

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).

**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

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 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.

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.

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.

**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.


As a final comment, since this manuscript is about FEM of bone scaffolds, the results of monotonic compression tests on bulk materials and scaffolds are, respectively, just an input to set the material properties in a finite element model and a comparison basis for the output for validation. Thus, since these are boundary conditions of the subject matter of the manuscript, required to reproduce FEM results, they are presented in this section.

#### *2.2. Methods*

Several numerical models have been developed in order to simulate the behaviour of the three cylindrical scaffolds mentioned in Section 3.1 when subjected to a monotonic compression test. To this objective, ANSYS v17 finite elements software has been employed.

The numerical models developed within the frame of this study follow the geometrical characteristics defined in Table 2 in terms of fibre diameter (0.25 mm), strand distance (1.25/1.00/0.75 mm) and layer height (0.20 mm). A diameter of 4 mm has also been considered for the samples. Regarding the total height of the scaffolds, due to geometrical considerations to enable the idealized interlayer link modelling, there is little difference between the real height (4.00 mm) and the height considered in the numerical models (4.05 mm).

Considering the double symmetry, not only of the geometry of the scaffolds, but also of the loadings to be considered, just one quarter of the scaffolds has been represented in the numerical models, as shown in Figure 7.

**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 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.

**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.

#### **3. Results**

#### *3.1. G25\_nanoHA 1.25*

In Figure 10a, the stress–strain curve for the G25\_nanoHA\_1.25 scaffold is shown. The strain (ε = Δl/l) is then calculated as the displacement of the upper plane divided by the total height of the scaffold, 4.05 mm. Additionally, the stress is calculated as the force to be applied in order to obtain a certain displacement of the upper plane divided by the cross-sectional area of the scaffold (π·ϕ2/4 = 4<sup>π</sup> mm2), where <sup>ϕ</sup> is the diameter of the cylindrical scaffold. In this regard, this stress can be considered as an "apparent stress". This is a very important difference, because the stress–strain relationship derived by this way will differ from the actual stresses obtained using the FE model, depicting stresses of every prismatic solid FE according to a colour scale.

**Figure 10.** Results of FEM simulation of a monotonic compression test on a G25\_nanoHA\_1.25 Scaffold: (**a**) Stress–strain curve; (**b**) Nodal stresses (Von Mises equivalent stress, in Pascals).

Thus, defining the elastic modulus of the scaffold E as the stress divided by the strain for a strain value of 0.03 (within this range, the stress–strain relationship is almost linear) leads to the following Equation (1):

$$E = \frac{0.1955}{0.03} = 6.60 \text{ MPa} \tag{1}$$

Additionally, in the following Figure 10b, the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown. Von Mises stress has been chosen for comparison because it is able to combine stresses in several axes in a single stress for comparison. Despite the fact there is a clearly dominant stress direction, since this is a simulated monotonic uniaxial compression test, it is combined with tangential and normal stresses in other directions at the spots where the fibres are crossing each other; therefore, simply plotting the principal stresses will underestimate the real triaxial stress state. As reflected in that figure, loads are mainly transferred through the columns formed by the intersections between perpendicular fibres.

It is noteworthy to remark the differences between the stresses of a stress–strain relationship (Figure 10a) and those depicted in the FEM (Figure 10b). The first going up to 0.5 MPa, while the second is going up to 13.2 MPa at substep 150. This is because, while the first is the quotient between the applied force divided by the cylinder area, the second is the actual stress at every SOLID 187 FE.

#### *3.2. G25\_nanoHA\_1.00*

Analogously, the same procedure is applied to the G25\_nanoHA\_1.00 numerical model. Accordingly, its stress–strain curve is shown in the following Figure 11a, with the corresponding Von Mises equivalent stress in Figure 11b. Moreover, the elastic modulus of the scaffold E, as the stress divided by the strain for a strain value of 0.03 leads to the following Equation (2):

$$E = \frac{0.2342}{0.03} = 7.90 \text{ MPa} \tag{2}$$

**Figure 11.** Results of FEM simulation of a monotonic compression test on a G25\_nanoHA\_1.00 Scaffold: (**a**) Stress–strain curve; (**b**) Nodal stresses (Von Mises equivalent stress, in 10−<sup>1</sup> Pascals).

It is noteworthy to remark the differences between the stresses in a stress–strain relationship (Figure 11a) and those depicted in the FEM (Figure 11b). The first going up to 0.28 MPa, while the second is going up to 20.4 MPa at substep 150. This is because, while the first is the quotient between the applied force divided by the cylinder area, the second is the actual stress at every SOLID 187 FE.

#### *3.3. G25\_nanoHA\_0.75*

Finally, the stress–strain curve of the G25\_nanoHA\_0.75 is shown in Figure 12a, while the corresponding equivalent Von Mises Stress is depicted in Figure 12b. The Young's modulus E is then calculated as the slope of the stress–strain curve at the point with 0.03 strain, as derived in the following Equation (3):

$$E = \frac{0.2342}{0.03} = 15.49 \text{ MPa} \tag{3}$$

**Figure 12.** Results of FEM simulation of a monotonic compression test on a G25\_nanoHA\_0.75 Scaffold: (**a**) Stress–strain curve; (**b**) Nodal stresses (Von Mises equivalent stress, in 10−<sup>1</sup> Pascals).

It is noteworthy to remark the differences between the stresses in the stress–strain relationship (Figure 12a) and those depicted in the FEM (Figure 12b). The first going up to 0.28 MPa, while the second is going up to 13 MPa at substep 150. This is because, while the first is the quotient between the applied force divided by the cylinder area, the second is the actual stress at every SOLID 187 FE.

#### **4. Discussion**

The results obtained using an FEM of the different scaffolds G25 nanoHA 1.25, 1.00 and 0.75 are compared to the corresponding monotonic compression tests on the Young's modulus basis (Table 3). As can be seen in all the cases, the finite element model with idealized geometry showed a Young's modulus higher than the one measured using direct testing. Moreover, this improvement was higher than the mean value plus 2–3 times the standard deviation of each scaffold. In the case of G25 nanoHA 1.25, the relationship of the FEM value to the mean measured value was 6.60/2.25 = 2.93 times higher, while in the case of the other two, G25 nanoHA 1.00 and G25 nanoHA 0.75, it was 1.54 and 2.45, respectively.

**Table 3.** Comparison between the Young's modulus obtained using FEM and direct scaffold testing, with corresponding Mean and Standard Deviation values.


Therefore, after repeating the FEM analysis several times, changing mesh sizes, contact definition and calculation steps and looking at the results, the difference was so high that it could not be explained by result scattering only. In fact, there are several assumptions on the FEM analysis on an idealized scaffold model that are not realistic at all, each contributing to these differences. For instance, a nominal fibre diameter, fibre position and straightness or strand distance are perfect in the model, but in the real world, with a 3D-printing process at such minimal scales, these parameters are impossible to be guaranteed and difficult to be controlled in real time. Thus, such imperfections are sources of uncertainty in terms of mechanical behaviour and the only way to deal with this is to set some execution tolerances and apply a design factor to the idealized geometry, as is common practice in structural engineering to deal with such execution imperfections.

Hence, the first step is to identify the main *prima facie* uncertainty sources. As such unintentional imperfections are impossible to reproduce, isolated from other errors and with accuracy in real bone scaffolds for testing, the best approach is to model and simulate such imperfections and see their influence on the results. In order to be able to simulate it, each uncertainty source is considered as an independent random variable, whose result is expected to be within a certain interval, with a reasonable grade of reliability, considering an according execution control procedure to guarantee some defined tolerances. Finally, all the variables need to be combined in a single design safety factor considering the global uncertainty, since each one will contribute to it to a certain degree.

As a last comment, a bone scaffold presents time-dependent properties in vivo. In fact, it undergoes deterioration processes related to fatigue, because a bone suffers from cyclic loadings and equivalent biological remodelling acting such as corrosion, since the scaffold material is progressively removed and substituted by bone tissue and this points for an initial deterioration up to a valley point where it starts to recover mechanical properties. Accordingly, the best approach is to develop a simplified model to take into account these two deterioration processes during service life, such as those presented by Calderon Uriszar-Aldaca et al. [50–52], developing an additional correction factor to take these into account.

#### *4.1. Sensitivity Analysis*

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

**Source Base Variation 1 Variation 2** Poisson's Ratio 0.4 0.375 0.35 Position of fibres At plane of sym. Out plane of sym. - Fibre diameter (mm) 0.25 0.235 0.22 Strand distance (mm) 1.25 1.2 1.3 Layer height (mm) 0.2 0.215 0.23 Sample diameter (mm) 4 3.9 4.1 Deflection of fibres (mm) 0 0.06 0.12 Straightness of columns Straight Alternated Cumulated Existence of broken fibres No broken fibres breakage at crossing far from crossing

**Table 4.** Comparison between base G25\_nanoHA\_1.25 characteristics and developed variations for sensitivity analysis.

#### 4.1.1. Poisson's Ratio of the Base Material

As mentioned in Section 3.1 materials, a Poisson's coefficient of 0.4 has been considered for the base material as common practice for polymeric materials. Nevertheless, as no specific test has been carried out in the frame of the project in order to determine this value accurately, it can be considered as a source of uncertainty.

Accordingly, three simulations have been carried varying the Poisson's ratio with 0.4, 0.375 and 0.35 as inputs for G25\_nanoHA\_1.25 scaffold finite element model. Then, the stress–strain relationship and the elastic modulus has been obtained following the same procedure to see how this variation influences the results.

Thus, Figure 13 shows the stress–strain relationship when varying the Poisson's ratio and the elastic modulus obtained as the slope of such stress–strain relationship, the curves are so close that they stack over each other. Moreover, Table 5 summarize the results. According to these results, the influence is almost negligible, i.e., after increasing the Poisson's ratio of the base material, only a slight increase in the elastic modulus of the scaffold can be appreciated, lowering only up to 99.09% in stiffness.

**Figure 13.** Influence of the Poisson's ratio on the stress–strain relationships from FEM.

**Table 5.** Compressive elastic modulus of a G25\_nanoHA\_1.25 scaffold for different values of the Poisson's ratio.


#### 4.1.2. Position of the Fibres within the Scaffold

According to Section 3 on results, the compression loads are mainly transferred to the basement through the columns formed by the intersections between the perpendicular fibres of the scaffold (Figure 10). Hence, due to the relatively low value of the sample diameter in comparison with the strand distance, the relative position of the fibres within the scaffold will affect the number and location of these columns, having some influence on the compressive stiffness of the scaffold, which can happen as a result of punching position variability.

As shown in Figure 14, two different configurations have been analysed. In the first one, at each layer of the scaffold, there is a fibre located in the relevant plane of symmetry. On the other hand, in the second configuration, at each layer of the scaffold, the relevant plane of symmetry is in the centre of the gap between two fibres.

**Figure 14.** Configurations considered to analyse the influence of the position of the fibres on the compressive stiffness of the scaffold (G25\_nanoHA 1.25).

Therefore, a detail of the geometry of the numerical models used for both configurations is shown in the following Figure 15:

**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).

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.


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.

**Figure 17.** Nodal stresses of G25\_nanoHA\_1.25 (Von Mises equivalent stress): (**a**) Configuration 1, in Pascals; (**b**) Configuration 2, in 10−<sup>1</sup> 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.


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).

**Figure 19.** Nodal stresses of G25\_nanoHA\_1.25 (Von Mises equivalent stress): (**a**) ϕfibre = 0.235 mm, in 10−<sup>1</sup> Pascals; (**b**) ϕfibre = 0.220 mm, in 10−<sup>1</sup> 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.

**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.


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).

**Figure 21.** Nodal stresses of G25\_nanoHA\_1.25 (Von Mises equivalent stress): (**a**) d = 1.20 mm, in 10−<sup>1</sup> Pascals; (**b**) d = 1.30 mm, in 10−<sup>1</sup> Pascals.

#### 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.


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 23.** Nodal stresses of G25\_nanoHA\_1.25 (Von Mises equivalent stress): (**a**) h = 0.215 mm, in 10−<sup>1</sup> Pascals; (**b**) h = 0.230 mm, in 10−<sup>1</sup> 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.

**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.


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 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.


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 27.** Nodal stresses of G25\_nanoHA\_1.25 (Von Mises equivalent stress): (**a**) f = 0.06 mm, in 10−<sup>1</sup> Pascals; (**b**) f = 0.12 mm, in 10−<sup>1</sup> 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:


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 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.

**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.

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.

**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).

**Figure 32.** *Cont*.

**Figure 32.** Influence of the straightness of the columns on the elastic modulus from FEM of G25\_nanoHA 1.25: (**a**) Influence on the stress–strain relationships of configuration 1; (**b**) Influence on the stress–strain relationships of configuration 2.

Accordingly, the values of the compressive elastic modulus of a G25\_nanoHA\_1.25 scaffold for each of the three different values considered in this analysis, for the deviation between the theoretical location of the fibres and their actual location, are shown in Table 12 for situation one and Table 13 for situation two.

**Table 12.** Compressive elastic modulus of a G25\_nanoHA\_1.25 scaffold for different values of the eccentricity or deviation in situation 1.


**Table 13.** Compressive elastic modulus of a G25\_nanoHA\_1.25 scaffold for different values of the eccentricity or deviation in situation 2.


As can be seen in Tables 12 and 13, the compressive elastic modulus of the scaffold with no deviation (6.70 MPa) differs from the compressive elastic modulus of the reference scaffold considered in the previous sections (6.60 MPa). This difference lies in the fact that, as mentioned, the numerical models used in both cases differ regarding the friction coefficient used for the contacts between the upper and lower planes and the scaffold fibres.

According to these results, when errors are cumulative (situation two), there is no significative influence, for the range of values considered, of the deviation between the theoretical location of the fibres and their actual location on the stress–strain relationship of a G25\_nanoHA\_1.25 scaffold subjected to a compression load. However, when errors alternate with respect to the theoretical location of the fibres (situation one), the influence of the deviation between the theoretical location of the fibres and their actual location on the stress–strain relationship of a G25\_nanoHA 1,25 scaffold subjected to a compression load may be considerable, because the straightness of the columns responsible for load transmission is much more affected in such situation.

#### 4.1.9. Existence of Broken Fibres

The G25\_nanoHA\_1.25 scaffold samples for testing have been printed using a theoretical fibre diameter of 0.25 mm, leading to relatively fragile fibres that can suffer damage during the printing process, the cooling and polymerization process or the tooling of the samples.

Thus, the influence of the existence of broken fibres on the stress–strain relationship of the G25\_nanoHA\_1.25 scaffold subjected to a compression load has been analysed using FEM, having considered the following alternatives:


All the fibre breakages considered in this study had a length of 0.15 mm and their location was randomly selected.

As in the previous section, in this case there is no symmetry in the geometry of the scaffolds; therefore, the whole scaffold has been represented in the numerical models. Therefore, in order to avoid lateral rigid body motions, which may impede the convergence of the numerical models, a friction coefficient of 0.3 has been considered in the contact between the upper plane and the upper fibres and the contact between the lower plane and the lower fibres. The results of this analysis are shown in Figure 33.

**Figure 33.** Influence of the existence of broken fibres on the stress–strain curve (G25\_nanoHA 1.25).

The value of the compressive elastic modulus of a G25\_nanoHA\_1.25 scaffold for each of the five different situations considered in this analysis concerning the existence of broken fibres is shown in Table 14.


**Table 14.** Compressive elastic modulus of a G25\_nanoHA\_1.25 scaffold for different situations concerning the existence of broken fibres.

As can be seen in Table 14, the compressive elastic modulus of the scaffold with no broken fibres (6.70 MPa) differs from the compressive elastic modulus of the reference scaffold considered in the previous sections (6.60 MPa). This difference lies in the fact that, as mentioned, the numerical models used in both cases differ regarding the friction coefficient used for the contacts between the upper and lower planes and the scaffold fibres.

According to these results, when the breakage of the fibres occurs far from any intersection between perpendicular fibres, there is no significant influence on the stress–strain relationship of a G25\_nanoHA 1.25 scaffold subjected to a compression load, as the paths for the load transmission (columns formed by the intersections between perpendicular fibres) are not affected.

On the other hand, some influence on the stress–strain relationship of a G25\_nanoHA\_1.25 scaffold subjected to a compression load can be seen when the breakage of the fibres is located in the intersection between perpendicular fibres. In this case, it is expected that the longer the breakage is and the higher the number of affected fibres is, the bigger the influence is.

#### *4.2. Design Safety Factor to Consider Uncertainty*

In view of the Section 4.1 outcome, it is clear that the great differences in terms of mechanical behaviour summarized in Table 3 could be explained by the printing process uncertainty. Any model will be performed according to an idealized geometry, conditions and material properties but, in the end, assuming this idealized behaviour is unsafe because 3D printing bone scaffolds at such small-scale leads to several manufacturing mistakes, which become uncertainty sources.

Nevertheless, although it is not possible to forecast exact predictions on mechanical behaviour while unable to perform a perfect 3D-printed bone scaffold, it is possible to perform enough safe ones. This procedure is analogous to that used in purely structural or mechanical engineering for other non-medical applications, where the forecasts are performed running idealized models with reduction factors to consider execution mistakes, material properties uncertainties, lack of homogeneity, etc.

Thus, such reduction factors need to be balanced to be conservative enough but not too penalizing for the structure and the building process. Therefore, several tolerances are defined for each execution variable that need to be checked and controlled during structure execution or, analogously in this case, for bone scaffold 3D printing. Hence, each variable is defined by a mean-measured value, μ, and a standard deviation, σ, assuming a Gaussian distribution according to central limit theorem, since each variable is dependent on several uncontrolled factors, and defining the check and control intensity in a way that there is a small probability, typically 5%, of being under a lower boundary limit or exceeding an upper bound with a reliability level typically of 75%. As an additional comment, these fifth percentile values are known as characteristic values.

Accordingly, each uncertainty source is defined as a random variable x1, x2, ... , x9 and, since the influence of each isolated random variable on the elastic modulus is linear dependent and biunivocal. the value of each random variable corresponds with an uncertainty factor, Fi, that is defined as the quotient between the actual elastic modulus and the theoretical one if 3D printed perfectly. For instance, the fibre diameter is the variable x3, the nozzle diameter is 0.25 mm, and it theoretically prints at 0.25 mm but, depending

on the printing speed, the fibre diameter could be lower. Thus, if the 3D printing process was perfect, then the fibre diameter would be 0.25 mm and the elastic modulus would be E = 6.6 MPa with a corresponding uncertainty factor F3 = 6.6/6.6 = 1.0. However, if the printing speed was higher, then the fibre diameter would be lower instead and, in the case that it was even 0.22 mm, then the elastic modulus would be E = 3.39 MPa and the corresponding factor F3 = 3.39/6.6 = 0.5136. Thus, because of the linear dependency, F3 is a random variable with a Gaussian distribution matching x3. Then, if we define a diameter tolerance and perform an execution control to guarantee that the fibre diameter is kept between 0.22 and 0.25 mm, with such an intensity that only 5% of the bone scaffolds present fibres below 0.22 mm and only 5% over 0.25 mm. Analogously, a mean, μi, and a standard deviation, σi, can be defined for each isolated uncertainty factor Fi, considering that in any Gaussian distribution, 90% of the probability is within the mean minus 1.6449 times the standard deviation and plus it. Hence, Table 15 summarizes the mean, μi, and standard deviation, σi, of each isolated uncertainty factor Fi, with corresponding limit values, Fi,min and Fi,max, for the interval.

**Table 15.** Uncertainty factors Fi and corresponding mean μ<sup>i</sup> and standard deviations σi, whose minimum and maximum values correspond to percentiles 5 and 95, respectively.


Thus, the next Figure 34 shows the probability density functions fi(xi) of each Gaussian random variable thus defined. For the sake of simplicity, each variable xi directly represents the uncertainty factor Fi = E(xi)/E(xtheoretical) defined as the quotient between the actual elastic modulus and the theoretical one if 3D printed perfectly.

**Figure 34.** Gaussian density function of each uncertainty reduction factor.

The mathematical expression of each Gaussian probability density function corresponds to Equation (4), as follows:

$$\mathbf{f\_i(x\_i)} = \frac{1}{\sigma\_{\mathbf{i}} \cdot \sqrt{2\pi}} \cdot \mathbf{e^{-\frac{(x\_i - \mu\_{\mathbf{i}})^2}{2 \cdot \sigma\_{\mathbf{i}}^2}}} \tag{4}$$

Hence, considering each isolated uncertainty variable as independent from each other, the probability density function of the multiplication is the multiplication of the independent probability density functions, see Equation (5) as follows:

$$\mathbf{f}(\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \cdots \cdot \mathbf{x}\_\bullet) = \mathbf{f}\_1(\mathbf{x}\_1) \cdot \mathbf{f}\_2(\mathbf{x}\_2) \cdot \cdots \cdot \mathbf{f}(\mathbf{x}\_\bullet) \tag{5}$$

Therefore, to derive the probability density function of the multiplication it is required a variable change. Accordingly, such variable change is shown in Equation (6), where *u* is the multiplication of eachN=9 variables *x*<sup>i</sup> and there are other N-1 variables zi.

$$\begin{cases} \mathbf{u} = \mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \cdots \cdot \mathbf{x}\_9\\ \mathbf{z}\_1 = \mathbf{x}\_1\\ \mathbf{z}\_2 = \mathbf{x}\_2\\ \vdots\\ \mathbf{z}\_8 = \mathbf{x}\_8 \end{cases} \tag{6}$$

Accordingly, the Jacobian determinant to execute the variable change is derived as the partial derivation of the variable multiplication by each variable, see Equation (7).

$$
\frac{\partial(\mathbf{x}\_{1}\cdot\mathbf{x}\_{2}\cdot\cdots\cdot\mathbf{x}\_{8})}{\partial(\mathbf{u},\mathbf{z}\_{1}\cdot\mathbf{z}\_{2}\cdot\cdots\cdot\mathbf{z}\_{8})} = \begin{vmatrix}
\frac{\partial\mathbf{u}\_{1}}{\partial\mathbf{u}} & \frac{\partial\mathbf{u}\_{1}}{\partial\mathbf{z}\_{1}} & \cdots & \frac{\partial\mathbf{u}\_{1}}{\partial\mathbf{z}\_{8}} \\
\frac{\partial\mathbf{u}\_{2}}{\partial\mathbf{u}} & \frac{\partial\mathbf{u}\_{2}}{\partial\mathbf{z}\_{1}} & \cdots & \frac{\partial\mathbf{u}\_{2}}{\partial\mathbf{z}\_{8}} \\
\vdots & \vdots & \vdots & \vdots \\
\frac{\partial\mathbf{u}\_{3}}{\partial\mathbf{u}} & \frac{\partial\mathbf{u}\_{3}}{\partial\mathbf{z}\_{1}} & \cdots & \frac{\partial\mathbf{u}\_{3}}{\partial\mathbf{z}\_{8}}
\end{vmatrix} = \begin{vmatrix}
0 & 1 & \cdots & \frac{\partial\mathbf{u}\_{1}}{\partial\mathbf{z}\_{8}} \\
0 & 0 & \cdots & \frac{\partial\mathbf{u}\_{2}}{\partial\mathbf{z}\_{8}} \\
\vdots & \vdots & \ddots & \vdots \\
\frac{1}{\mathbf{z}\_{1}\cdot\mathbf{z}\_{2}\cdots\mathbf{z}\_{8}} & \frac{-\mathbf{u}}{\mathbf{z}\_{1}\cdot\mathbf{z}\_{2}\cdots\mathbf{z}\_{8}} & \cdots & \frac{-\mathbf{u}}{\mathbf{z}\_{1}\cdot\mathbf{z}\_{2}\cdots\mathbf{z}\_{8}}
\end{vmatrix} = \frac{1}{\mathbf{z}\_{1}\cdot\mathbf{z}\_{2}\cdots\mathbf{z}\_{8}} \tag{7}
$$

Thus, the probability density function of the nine variables, once changed, is shown in Equation (9), as follows:

$$\mathbf{g}(\mathbf{u}, \mathbf{z}\_1, \cdots, \mathbf{z}\_8) = \mathbf{f}\_1(\mathbf{z}\_1) \cdot \mathbf{f}\_2(\mathbf{z}\_2) \cdot \cdots \cdot \mathbf{f}\_9(\frac{\mathbf{u}}{\mathbf{z}\_1 \cdot \mathbf{z}\_2 \cdot \cdots \cdot \mathbf{z}\_9}) \cdot \frac{1}{\mathbf{z}\_1 \cdot \mathbf{z}\_2 \cdot \cdots \cdot \mathbf{z}\_8} \tag{8}$$

In order to obtain the probability density function, the marginal distribution for variable *u*, which is obtained by the direct integration of the remaining variables z1 to z8, is required. See Equation (9) as follows:

$$\log(\mathbf{u}) = \iiint \cdots \int\_{\mathbb{R}^8} \left[ \mathbf{f}\_1(\mathbf{z}\_1) \cdot \mathbf{f}\_2(\mathbf{z}\_2) \cdots \mathbf{f}\_9 \left( \frac{\mathbf{u}}{\mathbf{z}\_1 \cdot \mathbf{z}\_2 \cdots \mathbf{v}\_8} \right) \cdot \frac{1}{\mathbf{z}\_1 \cdot \mathbf{z}\_2 \cdots \mathbf{v}\_8} \right] \cdot \mathbf{dz}\_1 \cdot \mathbf{dz}\_2 \cdots \cdot \mathbf{dz}\_8 \tag{9}$$

As a last comment on Equation (9), it is worthy to remember that each fi(xi) is a normal Gaussian distribution according to Equation (4) and that, despite being multiplied by the Jacobian fraction at the end, it is not directly integrable. Therefore, any solution can be obtained only by complex numerical means. Nevertheless, there is another simpler, more practical approach, considering that the distribution of the multiplication of Gaussian distributions is a Gaussian distribution itself. Then, it is possible to derive the mean and the standard deviation simply by taking into account their properties and their relationship with the expected value, E(x), and variance, Var(x).

#### 4.2.1. Derivation of the Global Mean

The expected value corresponds with the first order moment of the probability distribution. Therefore, the expected value of the multiplication of variables is the ninth

integration of the probability density function of the multiplication multiplied by each variable xi, between −∞ and ∞. See Equation (10), which is as follows:

$$\mathbb{E}[\hat{\mathbf{x}}\_1 \cdot \hat{\mathbf{x}}\_2 \cdot \cdots \cdot \hat{\mathbf{x}}\_9] = \iiint \cdots \cdot \int\_{\mathbb{R}^9} [\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \cdots \cdot \mathbf{x}\_9 \cdot \mathbf{f}\_{\mathbf{x}\_1 \cdot \hat{\mathbf{x}}\_2} \cdots \mathbb{x}\_9 (\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \cdots \cdot \mathbf{x}\_9) \cdot \mathrm{d}\mathbf{x}\_1 \cdot \mathrm{d}\mathbf{x}\_1 \cdots \cdots \cdot \mathrm{d}\mathbf{x}\_9] \tag{10}$$

Nevertheless, as the probability density function of the multiplication for a set of independent variables is the product of the isolated probability density function of each variable, see Equation (5), then the Equation (10) can be developed into Equation (11), which is as follows:

$$\mathbb{E}[\mathbf{\hat{x}}\_1 \cdot \mathbf{\hat{x}}\_2 \cdots \mathbf{\hat{x}}\_9] = \iiint \cdots \int\_{\mathbb{R}^9} [\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdots \mathbf{x}\_9 \cdot \mathbf{f}\_{\mathbb{R}\_1}(\mathbf{x}\_1) \cdot \mathbf{f}\_{\mathbb{R}\_2}(\mathbf{x}\_2) \cdots \mathbf{f}\_{\mathbb{R}\_9}(\mathbf{x}\_9) \cdot \mathbf{dx}\_1 \cdot \mathbf{dx}\_1 \cdots \cdots \mathbf{dx}\_9] \tag{11}$$

Accordingly, as every integral is independent from each other, then they can be regrouped as in the following Equation (12), presenting them as the product of integrals:

$$\mathbb{E}[\mathbf{\hat{x}}\_1 \cdot \mathbf{\hat{x}}\_2 \cdots \mathbf{\hat{x}}\_\theta] = \int\_{-\infty}^\infty \mathbf{x}\_1 \cdot \mathbf{f}\_{\mathbf{k}\_1}(\mathbf{x}\_1) \cdot \mathbf{dx}\_1 \cdot \int\_{-\infty}^\infty \mathbf{x}\_2 \cdot \mathbf{f}\_{\mathbf{k}\_2}(\mathbf{x}\_2) \cdot \mathbf{dx}\_2 \cdot \cdots \cdot \int\_{-\infty}^\infty \mathbf{x}\_\theta \cdot \mathbf{f}\_{\mathbf{k}\_\theta}(\mathbf{x}\_\theta) \cdot \mathbf{dx}\_\theta \tag{12}$$

Thus, remembering that the mean is precisely the first order moment of a probability distribution, see Equation (13), which is as follows:

$$\int\_{-\infty}^{\infty} \mathbf{x}\_{\mathbf{i}} \cdot \mathbf{f}\_{\mathbb{k}\_{\mathrm{i}}}(\mathbf{x}\_{\mathbf{i}}) \cdot d\mathbf{x}\_{\mathbf{i}} = \mu\_{\mathrm{i}} \tag{13}$$

Then, by direct substitution, the mean of any variable multiplication is the product of the isolated variable means, see Equation (14), which is as follows:

$$\mathbb{E}[\mathbb{A}\_1 \cdot \mathbb{A}\_2 \cdot \cdots \cdot \mathbb{A}\_9] = \mu\_1 \cdot \mu\_2 \cdot \mu\_3 \cdot \cdots \cdot \mu\_9 \tag{14}$$

#### 4.2.2. Derivation of the Global Standard Deviation

Regarding the variance of the variable multiplication, it corresponds to the expected value of the squared product minus the squared expected value of the product, see (15), which is as follows:

$$\text{VAR}(\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \mathbf{x}\_3 \cdot \dots \cdot \mathbf{x}\_9) = \text{E}\left[ (\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \mathbf{x}\_3 \cdot \dots \cdot \mathbf{x}\_9)^2 \right] - \text{E}\left( [\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \mathbf{x}\_3 \cdot \dots \cdot \mathbf{x}\_9] \right)^2 \tag{15}$$

Now, making use of the expected value properties, the expected value of the product of variables is the multiplication of the expected value of such variables. Therefore, substituting in Equation (15) there it comes Equation (16), which is as follows:

$$\text{VAR}(\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \mathbf{x}\_3 \cdot \dots \cdot \mathbf{x}\_9) = \text{E}\left[\mathbf{x}\_1^2\right] \cdot \text{E}\left[\mathbf{x}\_2^2\right] \cdot \dots \cdot \text{E}\left[\mathbf{x}\_9^2\right] - \left(\text{E}[\mathbf{x}\_1]\right)^2 \cdot \left(\text{E}[\mathbf{x}\_2]\right)^2 \cdot \dots \cdot \left(\text{E}[\mathbf{x}\_9]\right)^2 \tag{16}$$

Now, considering that the expected value of a variable is the mean value E[xi] = μi, the expected value of the squared variable corresponds to Equation (17), which is as follows:

$$\mathbb{E}\left[\mathbf{x}\_{\mathrm{i}}^{2}\right] = \left(\mu\_{\mathrm{i}}^{2} + \sigma\_{\mathrm{i}}^{2}\right) \tag{17}$$

Then, by a simple substitution in Equation (16), Equation (18) is then derived as follows:

$$\text{VAR}(\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \mathbf{x}\_3 \cdot \cdots \cdot \mathbf{x}\_9) = \left(\mu\_1^2 + \sigma\_1^2\right) \cdot \left(\mu\_2^2 + \sigma\_2^2\right) \cdot \cdots \cdot \left(\mu\_9^2 + \sigma\_9^2\right) - \mu\_1^2 \cdot \mu\_2^2 \cdot \cdots \cdot \mu\_9^2 \tag{18}$$

4.2.3. Derivation of the Safety Factor

Thus, according to previous Equations (14) and (18), when substituting the mean and standard deviation values for the uncertainty factors, taken as isolated variables, which are summarized Table 15, the mean and standard deviation of the product of variables is derived. See Equations (19) and (20), which are as follows:

$$\mu = \mathbb{E}[\hat{\mathbf{x}}\_1 \cdot \hat{\mathbf{x}}\_2 \cdots \cdot \hat{\mathbf{x}}\_9] \approx 0.4578 \tag{19}$$

$$
\sigma = \sqrt{\text{VAR}(\mathbf{x}\_1 \cdot \mathbf{x}\_2 \cdot \mathbf{x}\_3 \cdot \cdots \cdot \mathbf{x}\_9)} \approx 0.1346\tag{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.

**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:

$$\mathcal{E}\_{\rm d} = \frac{\mathcal{E}\_{\rm M}}{\gamma\_{\rm MP}} \tag{21}$$

$$\mathbf{E\_d} = 0.4578 \cdot \mathbf{E\_M} \approx \frac{\mathbf{E\_M}}{2.2} \tag{22}$$

$$\mathbf{E\_d} = 0.2364 \mathbf{\dot{\cdot}} \mathbf{E\_M} \approx \frac{\mathbf{E\_M}}{4.25} \tag{2.3}$$

Additionally, this design safety factor can be lowered if the execution control becomes more intense; thus 99% of the bone scaffolds presents the variables within the defined tolerance limits instead 95%, or if the intensity remains the same but the tolerance limits are lowered to ensure that geometry deviations are mitigated and derived, the uncertainty factors are lowered accordingly. The fibre diameter, layer height and straightness of the columns are key uncertainty sources affecting the bone scaffold behaviour. Therefore, the best strategy to improve the resulting scaffold behaviour is to improve the execution control and reduce the tolerances in these three variables.

#### 4.2.4. Use Case

For instance, let us have a use case analogous to that presented in Figure 1, where a hollow cylindric scaffold made of G25 nanoHA 1.25, G25 nanoHA 1.00 or G25 nanoHA 0.75 is substituting some damaged length of a femur for tissue regeneration, as representative of a uniaxial loading case. The specific patient body weight is around 75 kg; therefore, the scaffold section must face F = 735 N of loading peaks. If the cylindric scaffold length was 20 mm, the external diameter was 50 mm and the inner diameter was 10 mm, then the cross-sectional area would be A = 1885 mm2. Thus, the axial apparent stress of the scaffold would be σ<sup>x</sup> = F/A = −0.39 MPa, with σ<sup>y</sup> = σ<sup>z</sup> = 0 Mpa, since this is a monoaxial loading and a cylindric loading state, and the axial strain ε<sup>x</sup> is, according to Hooke's Equation (24) and taking into account the Poisson coefficients μyx and μzx, as follows:

$$
\varepsilon\_{\mathbf{x}} = \frac{\sigma\_{\mathbf{x}}}{\mathbf{E}\_{\mathbf{x}}} - \mu\_{\mathbf{yx}} \cdot \frac{\sigma\_{\mathbf{y}}}{\mathbf{E}\_{\mathbf{y}}} - \mu\_{\mathbf{zx}} \cdot \frac{\sigma\_{\mathbf{z}}}{\mathbf{E}\_{\mathbf{z}}} = \frac{\sigma\_{\mathbf{x}}}{\mathbf{E}\_{\mathbf{x}}} \tag{24}
$$

Now, using the apparent young modulus obtained from the idealized models EA,FEM already presented in Table 3 without any correction factor, the mean young modulus, Em,d, does not present any change. Thus, the mean strain ε<sup>m</sup> is obtained by direct application of Equation (24). Then, the displacement δ is obtained by multiplying the strain by the scaffold length. Moreover, the real young modulus, EFEM, taking into account the real stresses of solid 187 elements at the scaffold, multiplied by the strain ε<sup>m</sup> 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).


Now, if the design limits were σmax = 40 MPa and δmax = 1.5 mm, then the selected scaffold type could be G25 nanoHA 0.75 or G25 nanoHA 1.25, as they are the only ones fulfilling both requirements regarding the maximum displacement and stress. Nevertheless, if the scaffold was finally made of such materials, considering that the realistic scaffold will present lower stiffness due to several uncertain sources, then the calculation needs to be readjusted, taking into account the mean expectable behaviour of the scaffold (Table 17). This can be made by applying the reduction factor γMP = 2.2 to derive a mean Young's modulus for design Em,d and repeating the calculations with this adjusted elasticity. As can be seen, by calculating it this way, the only scaffold type left fulfilling both requirements is G25 nanoHA 0.75.


**Table 17.** Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with mean correction factor (γMP = 2.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 Ek,d and repeating the calculations with this adjusted elasticity. As can be seen, calculating it this way discards the three initial designs, since the three scaffolds present a stress over 40 MPa (Table 18).

**Table 18.** Derivation of displacement δ and stress σ at hollow cylindric scaffold within a femur bone with characteristic correction factor (γMP = 4.25).


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:


Reducing the strand distance from 0.75 to 0.5 leads to an increase in the amount of "resisting columns", from 21 to 45; this results in an expectable stress of 37.47 MPa, taking into account the previous results in columns made of spatially crossed cylinders of the same geometry. Moreover, increasing the fibre diameter instead, to pass from 60.2 to 40 MPa, will require a cross sectional augmentation by a scale factor equal to the square root of the quotient F = (60.2/40)0.5, giving a 0.307-mm diameter. Nevertheless, the diameter augmentation reduces the fibre cross-sectional curvature, increasing contact in a non-linear basis; therefore, it should be enough to increase the diameter up to 0.3 mm.

Thus, this patient will require a G25 nanoHA 0.75 with a fibre diameter of 0.3 mm and a G25 nanoHA 0.5 with fibre diameter of 0.25 mm, both produced from 55% PEOT/PBT + 45% nanoHA material or, finally, a G25 nanoHA 0.75 with 0.25 mm fibre diameter produced from an improved material.

#### **5. Conclusions**

Thus, comparing the results of finite element modelling to the bone scaffold compression tests according to considered methodology, and once the uncertainty sources derived from manufacturing process are analysed during discussion, the following conclusions can be summarized:

A methodology to perform FEM of 3D-printed bone scaffolds of any initial bulk material properties and geometry is defined. Moreover, it is experimentally calibrated and validated by direct testing of the scaffolds. The bone scaffold is then represented as an equivalent finite element to enable the FEM of bone tissue discretized with a compatible mesh, enabling the required characteristics to be studied, depending on bone topology.

The methodology is executed for a 55% PEOT/PBT + 45% nanoHA bulk material and three different scaffold geometries. The results coming from the FEM and the monotonic compression test of the bone scaffolds show significative differences, being that the idealized simulation is always better than the actual bone scaffolds. Considering the tests' mean values and standard deviations, the differences cannot be explained by result scattering only.

A deep look at the scaffolds show several fabrication imprecisions, meaning a resulting geometry far different from the theoretical one considered for FEM, such as different strand distance, layer height, fibre diameter, Poisson's ratio, fibre position, sample diameter, straightness of the columns, existence of broken fibres. This is due to the small imprecisions adding up layer-by-layer as large clinically relevant scaffolds (i.e., volume > 1 cm3) are being built. Thus, there are up to nine uncertainty sources detected in scaffolds coming from different variables.

Consequently, a sensitivity analysis is performed for each uncertainty source, deriving the influence or uncertainty factor and reducing the expected behaviour of each isolated source in terms of scaffold stiffness. Then, some tolerances are defined for each variable to be controlled during execution, with such an intensity that every variable presents a 95% probability of being within the selected interval in a Gaussian distribution.

The two main effects that can give strong changes to the elastic modulus are the strand overlapping and the fibre diameter. The first can be accommodated by design planning the shifting of the planes while printing, and the second may be modulating the printing speed, these will contribute to palliate the effect and minimize the corresponding design factor. A global design safety factor is defined in a predictive and conservative scenario, combining each isolated uncertainty factor in a single value taken from the mean value and fifth percentile of the Gaussian distribution of the product of each reducing uncertainty factor.

Thus, the aprioristic FEM of a bone tissue region with the inclusion of bone scaffold volume made of a certain bulk material becomes possible simply by using solid 185 equivalent finite elements to mesh it, assigning as material properties the ones derived from the FEM simulation of the monotonic tests of the scaffold, divided by the global design safety factor. This enables the further topological optimization of the scaffold along the bone tissue and the material selection, depending on the requirements.

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

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

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding authors.

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

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

