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

Abstract

Additive manufacturing (AM) of scaffolds enables the fabrication of customized patient-specific implants for tissue regeneration. Scaffold customization does not involve only the macroscale shape of the final implant, but also their microscopic pore geometry and material properties, which are dependent on optimizable topology. A good match between the experimental data of AM scaffolds and the models is obtained when there is just a few millimetres at least in one direction. Here, we describe a methodology to perform finite element modelling on AM scaffolds for bone tissue regeneration with clinically relevant dimensions (i.e., volume > 1 cm³). 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³, the simulations overestimated the scaffold performance, resulting in 150–290% stiffer than average values obtained in the validation tests. The results mismatch highlighted the relevance of the lack of printing accuracy that is characteristic of the additive manufacturing process. Accordingly, a sensitivity analysis was performed on nine detected uncertainty sources, studying their influence. After the definition of acceptable execution tolerances and reliability levels, a design factor was defined to calibrate the methodology under expectable and conservative scenarios.

1. Introduction

In order to manufacture patient-specific bone implants for tissue regeneration (Figure 1), it is important to keep in mind that any bone is a composite structure subjected to constant evolution normally composed of a high percentage of inorganic composition and a smaller organic fraction. The inorganic part is present intrafibrillarly in the bone and occupies up to 40% by volume filling with various degrees of space around the mineralized fibres, forming pore networks. Hence, human bones can be composed of high amounts of inorganics by volume, typically ranging from 50–60% up to 80–90% for some highly mineralized tissues, but always keeping some space for organics [1].

Consequently, with these boundary conditions, the engineering of patient-specific bone implants for tissue regeneration follows several approaches [2,3]. One of the most promising consists of the additive manufacturing of a temporary structure also known as a “scaffold”, suitable for cell growth, with required porosity to stimulate osteogenesis, as stated by Jakus et al. [4], Karageorgiou and Kaplan [5] and Hutmacher [6]. These temporary scaffolds must meet several requirements, such as biocompatibility, biodegradability and suitable mechanical properties [2,3,7,8,9,10,11,12,13], but also printability [14,15]. The mechanical properties are of great importance, since the scaffolds should ideally have properties matching those of the surrounding tissue, ensuring its suitability to bear mechanical loading similar to the original tissue. The mechanical properties are mainly dictated by the bulk material properties. However, the geometry and porosity of the fabricated scaffolds also play an important role in the final mechanical properties of the implant [16,17].

There are diverse materials being investigated, spanning from synthetic and natural polymers, ceramics and combinations of these, with the research focused on finding an optimal formulation for a successful stimulation of bone healing [2,3], enabling the improvement of implants suitable for each patient condition in the future. Thereby, additive manufacturing of a bone scaffold requires a compatible biocomposite with reinforcing inorganic fillers to mimic natural bone composition and structure, while improving its mechanical behaviour to resist service life loadings [18,19]. Accordingly, there are many studies regarding material alternatives [20], such as silicate nanocomposites [21,22], graphene nanocomposites [23] and the extensive work on hydroxyapatite or nano-apatite [24,25,26,27,28]. Additionally, other studies focused on the mechanical optimization and performance of such material alternatives in the presence of micro and nano hydroxyapatite [29], or on bone regeneration and thermomechanical properties [30,31].

In this manuscript, a generalized procedure is developed for the numerical modelling of additive manufactured scaffolds for bone regeneration with clinically relevant dimensions (i.e., volume > 1 cm³). The procedure is applicable to scaffolds made of any previously mechanically characterized material, and this procedure is then specifically validated with poly(ethylene oxide terephthalate)/poly(butylene terephthalate) (PEOT/PBT) with nano-hydroxyapatite fillers. This generalized procedure allows the scaffolds to be designed by function, fixing at first the requirements of geometry and the stresses to be supported and then selecting the optimal filler amount and porosity. A stiffer scaffold, for example, would require a greater amount of filler and less pores for stress to reduce the concentration spots. Consequently, this leads to design gradients and patterns for a single bone scaffold, which can be achieved using established topology optimization techniques.

Figure 1. 3D printed bone scaffolds for tissue regeneration: (a) long bone defect-shaped scaffold printed with PEOT/PBT alone; (b) Scaffolds presenting longitudinal gradient patterns of plasma polymerized film deposition, visualized using methylene blue staining [32].

Figure 1. 3D printed bone scaffolds for tissue regeneration: (a) long bone defect-shaped scaffold printed with PEOT/PBT alone; (b) Scaffolds presenting longitudinal gradient patterns of plasma polymerized film deposition, visualized using methylene blue staining [32].

The mechanical properties of the bone scaffolds are a key aspect for the development of PEOT/PBT composites and a fast-complete characterization of these bulk materials can be carried out using monotonic tensile and compression tests according to standards. Nevertheless, the equivalent testing of non-standardized bone scaffolds made of such materials is not as easy, for the scaffolds need to be previously manufactured with different geometries and architectures, which slow down their rapid characterization. Therefore, a modelling procedure suitable to predict the mechanical performance of bone scaffolds with a certain pore architecture and bulk material mechanical properties has been developed.

Finally, a set of scaffolds has been produced and mechanically characterised as a way to confirm the suitability of the developed modelling protocol and its validation. Nevertheless, the 3D printing of clinically relevant scaffolds with the required micro-scale precision is a challenging task and the manufacturing process generally produces unintended imperfections causing uncertainty because of result scattering. Hence, a sensitivity analysis and statistical study have been conducted in order to derive some design factors taking into account these issues.

The final output is a finite element modelling (FEM) procedure for 3D-printed bone scaffolds able to forecast their mechanical behaviour depending on the bulk material properties. This will enable the optimization of the material development process and perform topology optimization of such scaffolds for bone tissues. This is not a topology optimization necessarily based on evolutionary structural optimization (ESO) algorithms, removing material where it is not working at certain stress level; nor it is based on additive evolutionary structural optimization (AESO) algorithms, simply adding materials where the stresses are above certain stress level, since the global shape of the bone implant needs to be kept “custom made”, performing an external geometry for a specific patient. It is a topology optimization based on the scaffold properties’ distribution within such a patient custom made geometry, i.e., manufacturing implants with stronger bone scaffolds at stress concentration “hot spots” by varying layer height, strand distance, fibre diameter (changing extrusion nozzle or printing speed), bulk material properties, plasma treatments, etc.

There are several studies regarding the FEM of 3D bone scaffolds for tissue regeneration. For instance, there are several bone scaffold aspects being currently analysed in terms of scaffold design [33], mechanical behaviour during failure [34], or covering the bone tissue regeneration and the progressive change in properties [35,36]. Nevertheless, the scope of such studies is focused on the mechanical behaviour of the 3D bone scaffold itself. Yet, there is little research on the holistic methodology connecting the bulk material properties with the scaffold mechanical behaviour within a 3D-printed bone, with mechanical properties depending on topology and size effect [37,38,39]. Moreover, this kind of predictive tool is important to forecast the expected mechanical behaviour of the scaffold on-site, using an optimizable geometry or bulk material.

In order to do it with FEM, there are two basic steps required. The first one is to represent every scaffold as an equivalent finite element, characterized by an approximately cubic shape with eight nodes, each being able to undergo a displacement in three degrees of freedom. The second step is to discretize the bone geometry in a compatible mesh, made of such solid finite elements, whose material properties are defined to match the mechanical behaviour of the scaffold (Figure 2A).

Accordingly, a flowchart of the required works is presented in Figure 2B. It starts with the scaffold design at iteration i and the material tests for mechanical characterization under tension and compression. Then, an FEM of the bone scaffold is completed, and a monotonic compression test is simulated in that model. Thus, the results from the simulation can be compared to an actual scaffold compression test to calibrate and validate the FE model. This last step is required because the FEM is based on a theoretical or idealized geometry, whereas the real scaffolds manufactured at such scales present faults and uncertainties that need to be taken into account.

For the characterization of additive manufactured scaffolds, an FEM and support methodology has been carried out in this study, whose main objective is to be able to transfer the global mechanical properties of the scaffolds to solid eight-node finite elements, which will be used in FEM analyses. This will enable mechanical calculations to be addressed at a macroscopic level (complete bones, etc.), minimizing the computational costs, because it will not be necessary to take into account, in these models, the exact configuration of the scaffolds, but their equivalent mechanical properties. Additionally, this FEM of additive manufactured scaffolds has different properties, as equivalent finite elements can be further used to model mechanical property gradients simply by progressively changing such finite element material properties accordingly along a certain direction.

2. Materials and Methods

2.1. Materials

The bulk materials composing the scaffolds for bone tissue engineering were selected from amongst recently investigated new formulations based on combinations of a well characterized synthetic block copolymer of poly(ethylene oxide terephthalate)/poly(butylene terephthalate) (PEOT/PBT) with a diverse filler concentration [32,40,41]. For the improvement of the mechanical properties, the scaffolds produced from PEOT/PBT, including 45% of nano hydroxyapatite (nanoHA), were chosen due to their superior properties (Figure 3) and, therefore, were also used for the FE modelling. The mean mechanical behaviour obtained from the stress–strain curves was used (Figure 4). The stress–strain curves were obtained from the mechanical testing of the bulk materials according to the ASTM F2027-08 standard guide [42], prescribing ASTM D638-08 and ASTM D695-10 for tensile and compressive properties, respectively [43,44], and according also to international standard ISO 527 [45,46,47] (Figure 4). A Poisson’s ratio of 0.4 was taken as usual practice for polymers, although there is some evidence that it could be higher in the case of PEOT/PBT [16]. The effect of the bulk material’s Poisson’s ratio on scaffold mechanical properties was studied in a sensitivity analysis.

For the scaffolds, monotonic compression tests, according to ASTM F2150-02 and ISO 604, were carried out on cylindrical scaffolds [48,49] to study the effect of several parameters, namely strand distance, the diameter of the extrusion needle (which determines the scaffold fibre diameter), layer height and the structure of the scaffold (external geometry, fibre position, scaffold diameter, scaffold height). Figure 5).

Three variations of the scaffold were produced and tested, with three specimens each. The variations tested involved a decrease in the strand distance, meaning a denser, less porous scaffold, with a smaller inner volume of void spaces. The geometry of the different scaffold variations is summarized in

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

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

, with corresponding Mean and Standard Deviation (SD) values from the tests of three replicates each.

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

2.2. Methods

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

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

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

For the modelling of the fibres that conform to the scaffold, only one type of element (SOLID187) has been used. The SOLID187 element is a higher order 3D element with a quadratic displacement behaviour that is well suited to modelling irregular meshes. This element is defined by 10 nodes (tetrahedral element with nodes at vertices and mid edges) with three degrees of freedom at each node (translations in the nodal X, Y and Z directions). Moreover, for the modelling of the contacts, two additional types of elements have been automatically introduced by ANSYS, i.e., CONTA174 and TARGE170. CONTA174 is used to represent contact and sliding between 3D “target” surfaces and a deformable surface defined by this element, whereas the TARGE170 element is used to represent various 3D “target” surfaces for the associated CONTA174 contact elements. Hence, for demonstrative purposes, a detail of the meshing for G25_nanoHA_1.25 FEM with such finite elements is shown in Figure 8.

Regarding the boundary conditions, as the fibre diameter (0.25 mm) is larger than the layer height (0.20 mm), there will be stacking between perpendicular fibres. In these numerical models, it has been considered that there is a perfect bonding between fibres, as shown in Figure 9a. Additionally, a lower horizontal plane has been included in these numerical models to limit, by means of a contact, the vertical displacement of the lower fibres. This plane presenting infinite stiffness and its displacements are totally constrained, see Figure 9b. Finally, an upper horizontal plane has also been included in these numerical models in order to impose, by means of a contact, vertical displacement in the upper fibres. This plane presents infinite stiffness too and its lateral displacements are totally constrained (Figure 9c).

Finally, regarding the actions and type of analysis, in order to simulate a monotonic compression test, a vertical displacement has been imposed on the upper plane of the numerical models. This vertical displacement will be transferred to the scaffold by means of the contact defined between the upper plane and the upper fibres. Due to the negligible influence on the results, the self-weight of the scaffolds has not been considered in this study. Additionally, the need for evaluating the non-linearities related to both the material and the existence of contacts, has demanded the resolution of the numerical models developed in this study by means of a non-linear static analysis. In this resolution, the hypothesis of small displacements has also been taken into account.

3. Results

3.1. G25_nanoHA 1.25

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

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

$$E=\frac{0.1955}{0.03}=6.60 \mathrm{MPa}$$

(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 \mathrm{MPa}$$

(2)

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 \mathrm{MPa}$$

(3)

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. Comparison between the Young’s modulus obtained using FEM and direct scaffold testing, with corresponding Mean and Standard Deviation values.

Scaffold		G25 nanoHA 1.25			G25 nanoHA 1.00			G25 nanoHA 0.75
		Model	Mean	SD	Model	Mean	SD	Model	Mean	SD
Young’s Modulus	(MPa)	6.60	2.25	0.36	7.90	5.13	1.01	15.49	6.33	0.53

). As can be seen in all the cases, the finite element model with idealized geometry showed a Young’s modulus higher than the one measured using direct testing. Moreover, this improvement was higher than the mean value plus 2–3 times the standard deviation of each scaffold. In the case of G25 nanoHA 1.25, the relationship of the FEM value to the mean measured value was 6.60/2.25 = 2.93 times higher, while in the case of the other two, G25 nanoHA 1.00 and G25 nanoHA 0.75, it was 1.54 and 2.45, respectively.

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

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

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

4.1. Sensitivity Analysis

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

Table 4. Comparison between base G25_nanoHA_1.25 characteristics and developed variations for sensitivity analysis.

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

.

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. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the Poisson’s ratio.

Poisson’s Ratio (-)	Elastic Modulus (MPa)
0.400	6.60
0.375	6.57
0.350	6.54

summarize the results. According to these results, the influence is almost negligible, i.e., after increasing the Poisson’s ratio of the base material, only a slight increase in the elastic modulus of the scaffold can be appreciated, lowering only up to 99.09% in stiffness.

4.1.2. Position of the Fibres within the Scaffold

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

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

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

The influence of the relative position of the fibres within the scaffold on the stress–strain relationship of the G25_nanoHA 1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 16).

The values of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold of both configurations are shown in

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

(the compressive elastic modulus of the scaffold has been defined as the stress divided by the strain for a strain value of 0.03).

According to these results, there is some influence of the relative position of the fibres within the scaffold on the compressive elastic modulus of a G25_nanoHA_1.25 scaffold, with configuration one being stiffer than configuration two, only reaching 94.85% of the stiffness. Finally, the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 17 for both the configurations.

4.1.3. Fibre Diameter

As already disclosed in Table 1, scaffolds have been printed using a theoretical fibre diameter of 0.25 mm. Nevertheless, the accuracy of the printer in terms of deposition speed, positioning, temperature and material viscosity and variations in the material flow rate can lead to the real value of the fibre diameter being different to the theoretical one.

The influence of the fibre diameter on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed by means of FEM (Figure 18).

The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the fibre diameter is shown in

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

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

.

According to these results, when the variation of the strand distance is low, their influence on the stress–strain relationship of a G25_nanoHA_1.25 scaffold subjected to a compression load is almost negligible, i.e., increasing the strand distance causes only a slight decrease in the elastic modulus of the scaffold, up to 99.55% of the original stiffness.

As the compression loads applied on the scaffold will be mainly transferred through the columns formed by the intersections between the perpendicular fibres of the scaffold, variations in the strand distance will have a significant impact on the compressive elastic modulus of the scaffold, only when these variations lead to an increase/decrease in the number of those intersections.

For two different values of the strand distance (d = 1.20 mm and d = 1.30 mm), the Von Mises equivalent stresses for a scaffold strain of 0.03 are shown in Figure 21, for d = 1.25 mm (see also Figure 10).

4.1.5. Layer Height

The tested G25_nanoHA_1.25 scaffolds have been printed using a theoretical layer height of 0.20 mm. The influence of the layer height on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 22).

The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the strand distance is shown in

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

4.1.6. Sample Diameter

In the frame of this study, monotonic compression tests and dynamic compression tests have been carried out on cylindrical samples with a theoretical diameter of 4 mm, extracted from G25_nanoHA_1.25 scaffolds. Nevertheless, the accuracy of the tools used to extract the samples from the printed scaffolds can lead to a real value for the sample diameter different from the theoretical one. The influence of the sample diameter on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 24).

The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the sample diameter is shown in

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

The influence of the deflection of the printed fibres on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (in Figure 26).

The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis for the deflection of the printed fibres is shown in

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

4.1.8. Straightness of the Columns Responsible for Load Transmission

G25_nanoHA_1.25 scaffolds have been printed using a theoretical strand distance of 1.25 mm. However, the accuracy of the positioning system of the printer can lead to a real value for the strand distance different from the theoretical one.

Hence, the influence of the strand distance on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load was analysed in Section 4.1.4 on strand distance. In that case, it was supposed that there was some variation in the strand distance, but it was also considered that this variation remained constant within the whole scaffold. Under these circumstances, the columns formed by the intersections between the perpendicular fibres, which are responsible for the transmission of the compression load applied on the scaffold, remained vertically straight. Nevertheless, the x and y positioning in any printer present deviations; therefore, the real positions are x ± Δx and y ± Δy, meaning that some executions could present alternate positioning or show a displacement tendency.

In this section, on the other hand, the influence of the strand distance on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has also been analysed. In this case, this variation in the strand distance affected the straightness of those columns.

Thus, the following two situations have been analysed:

• Alternate variations in the strand distance (situation one);
• Cumulative variations in the strand distance (situation two).

In the first situation, there is a certain eccentricity e between the theoretical location of the fibres and their actual location. In addition, it has been considered that the fibre corresponding to a certain layer is located on the opposite side, regarding the theoretical location of the fibres, of that of the fibres corresponding to the preceding and the following layers, see Figure 28.

A detail of the geometry of the numerical model used to simulate this situation is shown in Figure 29. In this situation, the following two values have been considered for the deviation of the fibres: e = 0.025 mm and e = 0.05 mm.

In the second situation, there is also a certain eccentricity between the theoretical location of the fibres and their actual location. On the other hand, in this case, the deviation is cumulative, as shown in Figure 30.

A detail of the geometry of the numerical model used to simulate this situation is shown in Figure 31. In this situation, two values have also been considered for the deviation of the fibres, e = 0.005 mm and e = 0.009 mm.

In both situations, it has been considered that the deviation between the theoretical location of the fibres and their actual location affects the fibres oriented in both directions.

As can be seen in Figure 29 and Figure 31, in these cases, there is no any symmetry in the geometry of the scaffolds; therefore, the whole scaffold has been represented in the numerical models. Therefore, in order to avoid lateral rigid body motions, which may impede the convergence of the numerical models, a friction coefficient of 0.3 has been considered in the contact between the upper plane and the upper fibres and the contact between the lower plane and the lower fibres.

The influence of the deviation between the theoretical location of the fibres and their actual location, affecting the straightness of the columns formed by the intersections between the perpendicular fibres, on the stress–strain relationship of the G25_nanoHA_1.25 scaffold subjected to a compression load has been analysed using FEM (Figure 32).

Accordingly, the values of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the three different values considered in this analysis, for the deviation between the theoretical location of the fibres and their actual location, are shown in

Table 12. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the eccentricity or deviation in situation 1.

Deviation (mm)	Elastic Modulus (MPa)
0.000	6.70
0.025	5.76
0.050	3.99

for situation one and

Table 13. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different values of the eccentricity or deviation in situation 2.

Deviation (mm)	Elastic Modulus (MPa)
0.000	6.70
0.005	6.69
0.009	6.67

for situation two.

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

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

4.1.9. Existence of Broken Fibres

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

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

• No broken fibres within the scaffold—Situation one.
• One broken fibre within the scaffold (breakage located in one intersection between perpendicular fibres)—Situation two.
• Two broken fibres within the scaffold (breakages located in two different intersections between perpendicular fibres)—Situation three.
• One broken fibre within the scaffold (breakage located far from any intersection between perpendicular fibres)—Situation four.
• Two broken fibres within the scaffold (breakages located far from any intersection between perpendicular fibres)—Situation five.

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

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

The value of the compressive elastic modulus of a G25_nanoHA_1.25 scaffold for each of the five different situations considered in this analysis concerning the existence of broken fibres is shown in

Table 14. Compressive elastic modulus of a G25_nanoHA_1.25 scaffold for different situations concerning the existence of broken fibres.

Situation (#)	Elastic Modulus (MPa)
Situation 1	6.70
Situation 2	6.57
Situation 3	6.43
Situation 4	6.70
Situation 5	6.70

.

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 x₁, x₂, …, x₉ 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, F_i, 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 x₃, 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 F₃ = 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 F₃ = 3.39/6.6 = 0.5136. Thus, because of the linear dependency, F₃ is a random variable with a Gaussian distribution matching x₃. 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 F_i, 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. Uncertainty factors F_i and corresponding mean μ_i and standard deviations σ_i, whose minimum and maximum values correspond to percentiles 5 and 95, respectively.

Variable (xi)	Source (-)	Fi,min (#)	Fi,max (#)	μi (#)	σi (#)
x1	Poisson’s Ratio	0.9909	1.0000	0.9955	0.0028
x2	Position of the fibres	0.9485	1.0000	0.9742	0.0157
x3	Fibre Diameter	0.5136	1.0000	0.7568	0.1478
x4	Strand Distance	0.9955	1.0000	0.9977	0.0014
x5	Layer Height	0.6167	1.0000	0.8083	0.1165
x6	Sample Diameter	0.9545	1.0530	1.0038	0.0299
x7	Deflection of the fibres	0.9712	1.0000	0.9856	0.0088
x8	Straightness of the columns	0.5955	1.0000	0.7978	0.1230
x9	Existence of broken fibres	0.9597	1.0000	0.9799	0.0122

summarizes the mean, μi, and standard deviation, σi, of each isolated uncertainty factor F_i, with corresponding limit values, F_i,min and F_i,max, for the interval.

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

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

$${\mathrm{f}}_{\mathrm{i}}\left({\mathrm{x}}_{\mathrm{i}}\right)=\frac{1}{{\mathsf{\sigma }}_{\mathrm{i}}·\sqrt{2\mathsf{\pi }}}·{\mathrm{e}}^{-\frac{{\left({\mathrm{x}}_{\mathrm{i}}-{\mathsf{\mu }}_{\mathrm{i}}\right)}^2}{2·{\mathsf{\sigma }}_{\mathrm{i}}^2}}$$

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

$$\mathrm{f}\left({\mathrm{x}}_1·{\mathrm{x}}_2 \cdots {\mathrm{x}}_9\right)={\mathrm{f}}_1\left({\mathrm{x}}_1\right)·{\mathrm{f}}_2\left({\mathrm{x}}_2\right) \cdots {\mathrm{f}}_9\left({\mathrm{x}}_9\right)$$

(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 each N = 9 variables x_i and there are other N-1 variables z_i.

$$\{\begin{array}{l}\mathrm{u}={\mathrm{x}}_1·{\mathrm{x}}_2 \cdots {\mathrm{x}}_9\\ {\mathrm{z}}_1={\mathrm{x}}_1\\ {\mathrm{z}}_2={\mathrm{x}}_2\\ ⋮\\ {\mathrm{z}}_8={\mathrm{x}}_8\end{array}$$

(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 \left({\mathrm{x}}_1·{\mathrm{x}}_2 \cdots {\mathrm{x}}_9\right)}{\partial \left(\mathrm{u},{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8\right)}=\left|\begin{array}{cccc}\frac{\partial {\mathrm{x}}_1}{\partial \mathrm{u}}& \frac{\partial {\mathrm{x}}_1}{\partial {\mathrm{z}}_1}& \cdots & \frac{\partial {\mathrm{x}}_1}{\partial {\mathrm{z}}_8}\\ \frac{\partial {\mathrm{x}}_2}{\partial \mathrm{u}}& \frac{\partial {\mathrm{x}}_2}{\partial {\mathrm{z}}_1}& \cdots & \frac{\partial {\mathrm{x}}_2}{\partial {\mathrm{z}}_8}\\ ⋮& ⋮& ⋮& ⋮\\ \frac{\partial {\mathrm{x}}_9}{\partial \mathrm{u}}& \frac{\partial {\mathrm{x}}_9}{\partial {\mathrm{z}}_1}& \cdots & \frac{\partial {\mathrm{x}}_9}{\partial {\mathrm{z}}_8}\end{array}\right|=\left|\begin{array}{cccc}0& 1& \cdots & \frac{\partial {\mathrm{x}}_1}{\partial {\mathrm{z}}_8}\\ 0& 0& \cdots & \frac{\partial {\mathrm{x}}_2}{\partial {\mathrm{z}}_8}\\ ⋮& ⋮& \ddots & ⋮\\ \frac{1}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8}& \frac{-\mathrm{u}}{{\mathrm{z}}_1^2·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8}& \cdots & \frac{-\mathrm{u}}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8^2}\end{array}\right|=\frac{1}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8}$$

(7)

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

$$ℊ\left(\mathrm{u},{\mathrm{z}}_1,\cdots ,{\mathrm{z}}_8\right)={\mathrm{f}}_1\left({\mathrm{z}}_1\right)·{\mathrm{f}}_2\left({\mathrm{z}}_2\right) \cdots {\mathrm{f}}_9\left(\frac{\mathrm{u}}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_9}\right)·\frac{1}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8}$$

(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 z₁ to z₈, is required. See Equation (9) as follows:

$$ℊ\left(\mathrm{u}\right)={\displaystyle \iiint }\cdots {{\displaystyle \int }}_{{\mathrm{R}}^8}\left[{\mathrm{f}}_1\left({\mathrm{z}}_1\right)·{\mathrm{f}}_2\left({\mathrm{z}}_2\right) \cdots {\mathrm{f}}_9\left(\frac{\mathrm{u}}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_9}\right)·\frac{1}{{\mathrm{z}}_1·{\mathrm{z}}_2 \cdots {\mathrm{z}}_8}\right]·{\mathrm{dz}}_1·{\mathrm{dz}}_2 \cdots {\mathrm{dz}}_8$$

(9)

As a last comment on Equation (9), it is worthy to remember that each f_i(x_i) 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 x_i, between −∞ and ∞. See Equation (10), which is as follows:

$$\mathrm{E}\left[{\hat{\mathrm{x}}}_1·{\hat{\mathrm{x}}}_2 \cdots {\hat{\mathrm{x}}}_9\right]={\displaystyle \iiint }\cdots {{\displaystyle \int }}_{{\mathrm{R}}^9}^{}\left[{\mathrm{x}}_1·{\mathrm{x}}_2 \cdots {\mathrm{x}}_9·{\mathrm{f}}_{{\hat{\mathrm{x}}}_1·{\hat{\mathrm{x}}}_2 \cdots {\hat{\mathrm{x}}}_9}\left({\mathrm{x}}_1·{\mathrm{x}}_2 \cdots {\mathrm{x}}_9\right)·{\mathrm{dx}}_1·{\mathrm{dx}}_1 \cdots {\mathrm{dx}}_9\right]$$

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

$$\mathrm{E}\left[{\hat{\mathrm{x}}}_1·{\hat{\mathrm{x}}}_2\cdots {\hat{\mathrm{x}}}_9\right]={\displaystyle \iiint }\cdots {{\displaystyle \int }}_{{\mathrm{R}}^9}^{}\left[{\mathrm{x}}_1·{\mathrm{x}}_2\cdots {\mathrm{x}}_9·{\mathrm{f}}_{{\hat{\mathrm{x}}}_1}\left({\mathrm{x}}_1\right)·{\mathrm{f}}_{{\hat{\mathrm{x}}}_2}\left({\mathrm{x}}_2\right)\cdots {\mathrm{f}}_{{\hat{\mathrm{x}}}_9}\left({\mathrm{x}}_9\right)·{\mathrm{dx}}_1·{\mathrm{dx}}_1\cdots {\mathrm{dx}}_9\right]$$

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

$$\mathrm{E}\left[{\hat{\mathrm{x}}}_1·{\hat{\mathrm{x}}}_2\cdots {\hat{\mathrm{x}}}_9\right]={{\displaystyle \int }}_{-\infty }^{\infty }{\mathrm{x}}_1·{\mathrm{f}}_{{\hat{\mathrm{x}}}_1}\left({\mathrm{x}}_1\right)·{\mathrm{dx}}_1·{{\displaystyle \int }}_{-\infty }^{\infty }{\mathrm{x}}_2·{\mathrm{f}}_{{\hat{\mathrm{x}}}_2}\left({\mathrm{x}}_2\right)·{\mathrm{dx}}_2\cdots {{\displaystyle \int }}_{-\infty }^{\infty }{\mathrm{x}}_9·{\mathrm{f}}_{{\hat{\mathrm{x}}}_9}\left({\mathrm{x}}_9\right)·{\mathrm{dx}}_9$$

(12)

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

$${{\displaystyle \int }}_{-\infty }^{\infty }{\mathrm{x}}_{\mathrm{i}}·{\mathrm{f}}_{{\hat{\mathrm{x}}}_{\mathrm{i}}}\left({\mathrm{x}}_{\mathrm{i}}\right)·{\mathrm{dx}}_{\mathrm{i}}={\mathsf{\mu }}_{\mathrm{i}}$$

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

$$\mathrm{E}\left[{\hat{\mathrm{x}}}_1·{\hat{\mathrm{x}}}_2\cdots {\hat{\mathrm{x}}}_9\right]={\mathsf{\mu }}_1·{\mathsf{\mu }}_2·{\mathsf{\mu }}_3 \cdots {\mathsf{\mu }}_9$$

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

$$\mathrm{VAR}\left({\mathrm{x}}_1·{\mathrm{x}}_2·{\mathrm{x}}_3 \cdots {\mathrm{x}}_9\right)=\mathrm{E}\left[{\left({\mathrm{x}}_1·{\mathrm{x}}_2·{\mathrm{x}}_3\cdots {\mathrm{x}}_9\right)}^2\right]-\mathrm{E}{\left(\left[{\mathrm{x}}_1·{\mathrm{x}}_2·{\mathrm{x}}_3 \cdots {\mathrm{x}}_9\right]\right)}^2$$

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

$$\mathrm{VAR}\left({\mathrm{x}}_1·{\mathrm{x}}_2·{\mathrm{x}}_3 \cdots {\mathrm{x}}_9\right)=\mathrm{E}\left[{\mathrm{x}}_1^2\right]·\mathrm{E}\left[{\mathrm{x}}_2^2\right] \cdots \mathrm{E}\left[{\mathrm{x}}_9^2\right]-{\left(\mathrm{E}\left[{\mathrm{x}}_1\right]\right)}^2·{\left(\mathrm{E}\left[{\mathrm{x}}_2\right]\right)}^2 \cdots {\left(\mathrm{E}\left[{\mathrm{x}}_9\right]\right)}^2$$

(16)

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

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

(17)

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

$$\mathrm{VAR}\left({\mathrm{x}}_1·{\mathrm{x}}_2·{\mathrm{x}}_3 \cdots {\mathrm{x}}_9\right)=\left({\mathsf{\mu }}_1^2+{\mathsf{\sigma }}_1^2\right)·\left({\mathsf{\mu }}_2^2+{\mathsf{\sigma }}_2^2\right) \cdots \left({\mathsf{\mu }}_9^2+{\mathsf{\sigma }}_9^2\right)-{\mathsf{\mu }}_1^2·{\mathsf{\mu }}_2^2 \cdots {\mathsf{\mu }}_9^2$$

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

$$\mathsf{\mu }=\mathrm{E}\left[{\hat{\mathrm{x}}}_1·{\hat{\mathrm{x}}}_2\cdots {\hat{\mathrm{x}}}_9\right]\approx 0.4578$$

(19)

$$\mathsf{\sigma }=\sqrt{\mathrm{VAR}\left({\mathrm{x}}_1·{\mathrm{x}}_2·{\mathrm{x}}_3 \cdots {\mathrm{x}}_9\right)}\approx 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.

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:

$${\mathrm{E}}_{\mathrm{d}}=\frac{{\mathrm{E}}_{\mathrm{M}}}{{\mathsf{\gamma }}_{\mathrm{MP}}}$$

(21)

$${\mathrm{E}}_{\mathrm{d}}=0.4578·{\mathrm{E}}_{\mathrm{M}}\approx \frac{{\mathrm{E}}_{\mathrm{M}}}{2.2}$$

(22)

$${\mathrm{E}}_{\mathrm{d}}=0.2364·{\mathrm{E}}_{\mathrm{M}}\approx \frac{{\mathrm{E}}_{\mathrm{M}}}{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.

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

$${\mathsf{\epsilon }}_{\mathrm{x}}=\frac{{\mathsf{\sigma }}_{\mathrm{x}}}{{\mathrm{E}}_{\mathrm{x}}}-{\mathsf{\mu }}_{\mathrm{yx}}·\frac{{\mathsf{\sigma }}_{\mathrm{y}}}{{\mathrm{E}}_{\mathrm{y}}}-{\mathsf{\mu }}_{\mathrm{zx}}·\frac{{\mathsf{\sigma }}_{\mathrm{z}}}{{\mathrm{E}}_{\mathrm{z}}}=\frac{{\mathsf{\sigma }}_{\mathrm{x}}}{{\mathrm{E}}_{\mathrm{x}}}$$

(24)

Now, using the apparent young modulus obtained from the idealized models E_A,FEM already presented in Table 3 without 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. 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. 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

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

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

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