Micro computed tomography based ﬁnite element models for elastic and strength properties of 3D printed glass scaffolds

In this study, the mechanical properties of glass scaffolds manufactured by robocasting are investigated through micro computed tomography ( μ − CT ) based ﬁnite element modeling. The scaffolds are obtained by printing ﬁbers along two perpendicular directions on parallel layers with a 90 ◦ tilting between two adjacent layers. A parametric study is ﬁrst presented with the purpose to assess the effect of the major design parameters on the elastic and strength properties of the scaffold; the mechanical properties of the 3D printed scaffolds are eventually estimated by using the μ − CT data with the aim of assessing the effect of defects on the ﬁnal geometry which are intrinsic in the manufacturing process. The macroscopic elastic modulus and strength of the scaffold are determined by simulating a uniaxial compressive test along the direction which is perpendicular to the layers of the printed ﬁbers. An iterative approach has been used in order to determine the scaffold strength. A partial validation of the computational model has been obtained through comparison of the computed results with experimental values presented in [10] on a ceramic scaffold having the same geometry. All the results have been presented as non-dimensional values. The ﬁnite element analyses have shown which of the selected design parameters have the major effect on the stiffness and strength, being the porosity and ﬁber shifting between adjacent layers the most important ones. The analyses carried out on the basis of the μ − CT data have shown elastic modulus and strength which are consistent with that found on ideal geometry at similar macroscopic porosity.


Introduction
The tissue engineering (TE) approach to bone repair aims at mimicking the biological and biomechanical environment found in the bone in order to drive progenitor cells towards a defined differentiation and obtain fully functional tissue as replacement at anatomical sites damaged by injuries or Excutive Editor: Ji-Zeng Wang B Pasquale Vena diseases [1]. Among material candidates, metals, polymers and ceramics are found [2]. The use of selected materials and the versatile additive manufacturing (AM) techniques bring to scaffold designs able to meet mechanical and morphological requirements for biocompatible materials [1,[3][4][5]. In particular, optimal micro-scale architecture with high interconnected porosity with the purpose to ensure cell penetration, vascular ingrowth and nutrient diffusion can be achieved.
As a representative example, Entezari et al. [6] presented ceramic bone scaffold in Strontium-hardystonite-Gahnite (Sr-HT-Gahnite) manufactured through 3D printing technology. The printing process showed to be flexible enough to design scaffold with different geometrical features (fiber diameter and spacing) spanning a wide range of macroscopic strength.
In the present work we deal with bioactive glass scaffolds [7] obtained through the robocasting technique [8]. A match between the suitability of the glass ceramic materials in cre-ating a stable tissue-device interface [9] with the capability of the manufacturing technique to tune morphological and mechanical properties [8,10] can be achieved.
Producing glass ceramic scaffolds featuring a strength, toughness and stiffness suitable to bear the mechanical burden under physiological conditions is a crucial factor. Indeed, unlike the high toughness exhibited by lamellar bone which constitutes the trabecular structures, the ceramic constituent material exhibits a brittle response [11]; in addition to that, microstructural defects such as micro-voids or microcracks which are intrinsic to the manufacturing technique substantially affect the overall strength of the scaffold.
Micro computed tomography(μ−CT ) is a non destructive technique which produces high resolution three dimensional model of the scanned sample [12], allowing for the identification of final microstructural features of the scaffold, and to identify defects as obtained from the printing technique. In addition to micro-morphological assessment, μ − CT is also powerful tool to assess the effect of micro-morphology of the scaffold on its mechanical properties.
Strength prediction for glass ceramic foams on the basis of μ − CT scans and analytical modeling of the brittle fracture failure mechanisms [13] allowed at determining an estimate of the scaffold strength as a function of the foam porosity. In this approach, tomography scans are used with the only purpose to determine the foam porosity.
μ − CT based finite element modeling, originally introduced to assess mechanical properties of trabecular bone [14], allows at estimating mechanical properties of the bone scaffolds allowing for a sample specific geometrical models [15][16][17]. However, although the estimate of the macroscopic elastic property of bone scaffold through finite element modeling is relatively well established (see for example [16]), the estimate of the scaffold strength is still a crucial point, which is not well consolidated. The simulation of the fracture process in brittle scaffolds is challenging; Miranda et al. [18] predicted the fracture strength of glass ceramic scaffold obtained through robocasting by using traditional finite element approach to strength assessment. They limited the analyses to ideal geometrical models i.e. not based on tomography scans. The strength was predicted by means of the critical stress criterion. More recently, Entezari et al. [6] made use of the extended finite element method (XFEM) [19] which resulted in an overall satisfactory prediction of the scaffold strength and the relevant fracture pattern.
The aim of this paper is to predict the macroscopic elastic modulus and the overall strength of glass ceramic scaffolds obtained through the robocast technique. A twofold objective is here proposed: i) a sensitivity analysis of the macroscopic mechanical properties with respect to the major geometrical parameters of the scaffold, including relative density (porosity), fiber spacing and fiber misalignment; ii) assessing the effect of the real microstructural features as obtained through the printing process on the macroscopic mechanical property. This latter aim is achieved through μ − CT based finite element models allowing for crack propagation in the brittle material. Strength normalized to the intrinsic material strength will be presented.

Scaffold fabrication by robocasting
Bioactive glass scaffolds based on a multicomponent meltderived silicate glass (47. were fabricated by robocasting, as described in detail elsewhere [8]. Briefly, glass particles with size below 32 µm were added to a solution of Pluronic F-127 (Sigma-Aldrich) for achieving a printable ink. Specifically, 35 vol.% of glass was embedded into an optically-clear water-based solution containing 27.5 wt.% of Pluronic F-127, which was previously stirred overnight while being maintained in an ice bath due to its thermosensitive behaviour. The glass-containing ink was then homogenized for 1 min by using a vortex mixer (Vibrofix VF1 electronic, Ika-Werk) at 2500 rpm and cooled for 1 min into an ice bath. Five mixing-cooling cycles were carried out to obtain good dispersion of glass particles [10,20]. Homogenized ink was cast into the printing cartridge connected to the pressurecontrolled robocasting machine (Tabletop-3Dn, nScrypt Inc., Orlando, FL, USA) and left to stabilize for 1 h before printing. The printing tower could move along the z axis only (vertical position) and the plate under the extrusion nozzle could move on the horizontal xy plane. Polymeric tips with an inner diameter of 410 µm (Nordson EDF Optimum®S-moothFlow™) were used to extrude the ink. The printing accuracy was 10 µm along the x and y axes and about 5 µm along the z axis. Scaffolds were printed on acetate sheets (Colour Laser Printer & Copier OHP Film, Folex AG, Seewn, Switzerland) that exhibited flatness and good adhesion with the wet ink, from which the substrate could then be easily detached after drying [21]. The printing software provided by nScrypt (Machine Tools 3.0) allowed every single movement of the printing head (raster pattern) to be programmed by compiling a design script similar to a text file. The printing speed was 2 mm/s and the pressure used to extrude the ink filaments was in the range of 1.24-1.51 bar. Bioactive glass scaffolds were printed in the form of porous cuboids (length = width = 7.5 mm) with a grid-like structure comprising 20 layers, each of which was 90 • -tilted with respect to the adjacent ones (total height about 4.5 mm). Once the printing process was completed, each scaffold was left to dry for 48 h in air and, then, sequentially treated in an electrically-heated furnace at 200 • C/30 min, 400 • C/30 min and 500 • C/30 min to allow organic phase removal, followed by glass sintering

− CT scans
In order to acquire 3D volume reconstructions, scaffolds were scanned by X-ray micro-computed tomography (1μ − CT ; Phoenix Nanotom S machine, General Electric Measurement and Control) using source voltage and current of 110 kV and 110 µA, respectively. No X-ray filters were used. A translational motion compensation was performed in order to ensure a perfect matching of the 0 • and 360 • shadow images. During μ − CT analysis, scaffold magnification of 10X was adopted; major scanning parameters included voxel size of 5.00 µm, rotation step of 0.50 • and exposure time of 1.5 s. The μ−CT projections were used to reconstruct the scaffold volume in 3D by means of the Radon transform [22] implemented in the software datos-x-reconstruction provided by the manufacturer. VG Studio Max 2.0 (Volume Graphics) was used for visualization and exportation of image stacks.

Assessing the effect of 3D printing parameters on scaffold properties
Three different printing parameters have been selected with the purpose to quantitatively assess the effect of the geometry on the scaffold stiffness and strength. More specifically, the following geometric parameters are investigated: i) interfiber distance (dx = dy = β D) in the printing plane, ii) fiber shifting between two parallel printing planes (s = αdx), iii) fiber compenetration between two adjacent printing planes c = γ D. In which, α, β and γ are three dimensionless parameters; γ = d−D D , with d is the distance between the center-lines of fibers belonging to two adjacent printing planes and D is the fiber diameter. A sketch showing the geometrical interpretation of the three parameters is shown in Fig. 1.
A voxel type finite element mesh [14] has been created for each of the generated geometries with the purpose to simulate a μ − CT based finite element mesh in which one single hexahedron element having a cubic shape is placed at each pixel location of the image.
To this purpose, the images from μ − CT scans where transformed in a 8-bit grey scale stack of images before filtering and thresholding. A Gaussian filter with two pixel radius was applied to reduce the image noise, then the Otsu algorithm was used to binarize the image sequence. As a last step a connectivity analyses (purify) was run with the purpose to remove isolated floating small bodies. All image processing was carried out in ImageJ [23]. The parametric analyses presented in this section are performed on a finite element model having an idealized geometry; this latter being obtained through a MATLAB script which simulates a pixellike binary image as those obtained after binarization of CT scans. Each cylindrical printed fiber is discretized with 30 pixels along the diameter and assuming a fictitious space resolution of 6 µm/pixel, the fibers exhibit a 180 µm diameter. The open source multigrid finite element code PAROSOL [24] has been used for the core solution of the elastic finite element analyses.
The fracture propagation process was simulated through an algorithm which annihilates the element stiffness in selected elements of the mesh on the basis of a stress based criterion.
The incremental process is a sequence of elastic analyses in which a uniaxial compression is simulated. According to Fig. 1, the nodes belonging to the top surface are subjected to a fictitious displacement u 0 z =ū along direction z, while displacement along the x and y directions are not constrained. The nodes belonging to the bottom surface are subjected to a displacement constrain u 0 z = 0. The rest of the boundary nodes are not subjected to any kinematic constraints.
This set of boundary conditions simulates a macroscopic strain z =ū L z in which L z is the total length of the domain of analysis along loading direction z. In order to restrain rigid body motions in the x-y plane three additional boundary conditions are applied in which the displacement along x and y directions are constrained on a node of the base, and displacement along y direction is constrained in another node of the base aligned along the x direction with the other constrained node.
The computational procedure is essentially based on a sequence of elastic analyses, all run with a prescribed compression displacement on the top base having a reference valueū = u 0 (see block diagram in Fig. 2), these analyses provide reference stress σ (ū) and strain (ū) fields; at each run a given number of element is annihilated. The criterion used for element annihilation is a Drucker Prager criterion able to simulate a tension/compression mismatch of the strength property of ceramic-based materials. The Drucker Prager equivalent stress σ D P is defined as: in which S i j is the deviatoric component and p = σ 11 +σ 22 +σ 33 3 is the hydrostatic component of the Cauchy stress tensor. If uniaxial tensile σ T and compressive σ C strength are given, the parameters α D P andσ D P are: here σ C is a positive parameter, being the absolute value of the compression stress at failure. The Drucker-Prager equivalent stress σ D P is calculated in the centroid of each of the cubic elements. All elements of the mesh are sorted for increasing σ D P and the top N toll most stressed elements are selected, in which N toll is a given number which defines the total amount of elements to be annihilated for the current iteration. In the list of N toll elements to be annihilated, the stress in the less stressed element defines the maximum local stress for that iteration. Being the less stressed element in the list the i −toll element in the mesh one defines the ratio The displacement for that specific loading increment is then scaled by (see bullets in the force displacement plot reported in the inset of Fig. 2); stress and strain fields are computed at the current increment by scaling reference values of stresses and strains accordingly, see flow-chart in Fig. 2. Under these conditions, all elements in the mesh will be subjected to an equivalent stress σ D P ≤σ D P . As the properties of an equivalent homogenized material are sought, the macroscopic stress Σ j zz is calculated at each loading step j: in which V is the macroscopic domain of analysis which includes the solid material and the empty spaces. Owing to the linearity of the model, here the local stress σ zz is scaled of the same amount ξ i−toll as done for the displacements. As the stress in the empty space is zero, the above integral reduces to: in which V s is the volume occupied by the solid phase only. The stress are provided in the centroid of the element; therefore, the integral in 6 reduces to: in which n e is the total number of elements in the mesh, v e is the volume of the single element (i.e. the pixel size to the third power) and σ i zz is the stress in the i-th element of the mesh. The macroscopic elastic modulus (Ē) along direction z is given by the simple relationship: in which Σ 0 zz is the macroscopic stress in the first loafing step in which no elements have been annihilated yet, and u 0 is the reference displacement applied on the top surface. The macroscopic material strength is obtained as the maximum value of Σ i zz achieved during the loading process.

Assessment of − CT based scaffold elasticity and strength
The 3D microstructure of the printed scaffold was quantitatively investigated through μ−CT by using a high resolution as reported in the section 2.2. The resulting pixel size is 5 µm. The 16bit acquired images were treated by means of a custom made analysis code in MATLAB coding environment. All images in the stack were transformed into 8bit images with 255 grey levels. The images were then treated with a thresholding algorithm as described in Ref. [25] and a black and white stack of images is eventually obtained. Four different subdomains have been selected for finite element analyses having cubic shape with 300 pixles on each edge (macroscopic volume 1.5 mm × 1.5 mm × 1.5 mm ). In each subdomain, 8 fiber intersections have been included (see Fig. 3). A connectivity filter was applied in order to remove unconnected regions through the bwlabeln built-in function in Matlab (connectivity parameter = 6). The Cartesian mesh was obtained by converting each solid voxel into an 8-noded hexahedral element [14,26] with isotropic linear elastic material properties. A fictitious value of the Young The size of the finite element domain for all subregions identified in the scaffold is 300 pixel×300 pixel×300 pixel.
For macroscopic elasticity and macroscopic strength estimation, uniaxial compression tests were simulated through a displacement controlled finite element run. The nodes belonging to the top surface are subjected to a fictitious displacement u 0 z =ū along direction z, while displacement along the x and y directions are not constrained. The nodes belonging to the bottom surface are subjected to a displacement constrain u 0 z = 0. The rest of the boundary nodes are not subjected to any kinematic constraints.
Strength has been calculated as the maximum macroscopic stress achieved in the incremental process, in which N toll = 1.3 × 10 −4 n e .

Elasticity and strength of scaffolds with ideal 3D printed geometry
Macroscopic uniaxial elastic moduli are reported in a dimensionless form asĒ = E F E M /E s in which E F E M is the macroscopic uniaxial elastic modulus as obtained through the post-processing of the simulated uniaxial compressive test and E s is the intrinsic elastic modulus of the solid phase in the scaffold. Similarly, strength of the scaffold is also reported in a dimensionless form asσ c = Σ zz /σ t , in which Σ zz is the macroscopic uniaxial compression strength as obtained by the pos-tprocessing of the finite element simulation of the compression test and σ t is the tensile strength of the solid material under tension. Figure 4 shows the dimensionless macroscopic stiffness and strength for different values of the shifting geometric parameter α. The obtained result is expected as increasing the shifting between fibers a progressively increasing role of the bending mechanisms occurs. As macroscopic porosity of all models obtained through fiber shifting is approximately the same (i.e φ ≈ 0.44), the decrease in stiffness and strength is mainly motivated through a change in mechanisms involved, with progressively more bending of the fibers as shifting increases. Figure 5 shows the damage distribution for two models with two values of fiber shifting: the model on the right panel is characterized by a lower value of shifting in comparison to the model shown in the left panel. The damage distribution in the model on the left shows that material failure originates at the interface between two perpendicular fibers where shear stress is expected to be substantial; while, the model on the right shows that, for higher value of fiber shifting, the material failure is dominated by the tensile stress rising on the tensile side of fibers under bending. Figure 6 shows the dimensionless macroscopic stiffness and strength for different values of the geometrical parameter β (i.e. the inter-axis fiber distance). The finite element simulations shows how macroscopic stiffness and strength substantially decrease by increasing the inter-axis distance (i.e. increasing the value of β). The macroscopic porosity of the models obtained through increasing parameter β increases from φ = 0.54 to φ = 0.78; this indicates that the decrease of material properties is mainly due to the increase of the porosity of the material. Figure 7 shows the dimensionless macroscopic stiffness and strength for different values of the geometrical parameter γ (i.e. the fiber superposition on adjacent layers). The stiffness and strength increase with the fiber superposition (i.e. with increasing the value of the geometrical parameter γ ). The macroscopic porosity of models obtained by increasing the fiber superposition ranges between the values of φ = 0.34 and φ = 0.63; in this range, the mechanical properties (both stiffness and strength) increase of a factor four (from 0.1 to 0.4, approximatively). This substantial increase of the stiffness and strength is given by the decrease of the porosity and by the increase of the contact area between two perpendicular fibers.
A partial validation of the proposed computational model is obtained by means of a comparison with the results presented in Ref. [6]. In Ref. [6], Sr-HT-Gahnite scaffolds were fabricated through the robocasting technique. The geometrical parameters describing the printed geometry are the same as those used in the present paper. In order to have comparable results, the intrinsic material tensile strength was selected so to match the experimental strength obtained through the computational prediction for the model exhibiting the lowest porosity (φ = 0.5, Σ T = 140 MPa). In order to have a match on the macroscopic strength at φ = 0.5, the intrinsic material strength of 495 MPa was set.
Once the intrinsic material tensile strength was chosen, three additional models with porosity 0.56, 0.57 and 0.60 where simulated. Figure 8 shows the strength as obtained through the proposed finite element model for four different values of

Elasticity and strength of scaffolds with − CT based 3D printed geometry
Four different subdomains were selected for macroscopic elasticity and strength assessment. In order to save computational costs, each subdomain contains eight fiber intersections. The domains are cubic in shape with approximately 200 pixel per side. One single finite element per pixel was used to create the finite element mesh. Uniaxial compression was simulated along the direction perpendicular to the plane of printing. Figure 3 shows a representative volume for μ − CT based finite element simulation of the compression tests. Figure 9 (top and bottom) shows a typical distribution of the material damage upon compression. Damage originating at the contact regions between two perpendicular fibers is shown. In addition to that, damage developing in the free central region is found. Figure 10 shows a typical stress strain response as obtained through the finite element simulation. The stress-strain relationship is mostly linear the until damage spreading reduces substantially the macroscopic stiffness. In order to reduce computational costs, simulations were not conducted until final failure. Table 1 shows the results for macroscopic normalized strength and stiffness for the four subdomains obtained from the μ − CT images of the scaffold.
The porosity of the four models were similar, ranging from φ = 0.37 to φ = 0.45. The model showing the lowest porosity consistently resulted in higher strength and stiffness.
In order to have a complete overview of the results obtained on the models with idealized geometry and a comparison of the μ − CT based models, macroscopic strength and stiffness versus macroscopic porosity for all models are represented in the Fig. 11. Strength and stiffness decreases with the porosity for all models, the strength and stiffness of the μ − CT based models fall within the global trend delineated by the models with idealized geometry. A third order polynomial fitting has been found to be suitable to closely describe both elasticity and strength. Both data have been fit using the equation P = a(1 − φ) 3 in which P is the macroscopic dimensionless mechanical property: i.e. the macroscopic elastic modulus and macroscopic compressive strength. Fitting parameters are a = 1.57, R 2 = 0.93 and a = 1.50, R 2 = 0.92 for stiffness and strength, respectively. It is worth noting that, for the dimensionless elasticity and strength fitting polynomial, the parameter a should be equal to 1. The results found in this fitting analysis indicate that the fitting function should not be used outside of the present porosity range.

Discussion
In this study a computational model to predict the macroscopic elastic and strength properties of 3D printed scaffold has been presented. In particular, the study is divided in two parts: the first part presents a sensitivity analysis with respect to geometrical parameters defining the printed geometry of the scaffold; the second part presents the application of the computational method to real 3D printed glass scaffolds. Figures 5, 6, 7 shows how fiber shifting, fiber spacing and fiber interpenetration affect the elastic and strength properties. It is worth noting that, despite the number of pixel in the diameter is not very high (60 pixels/elements) the sensitivity analyses, which has been presented for comparison purposes, clearly showed that each geometric parameter had Despite a wide porosity span can be obtained by varying different geometrical parameters (fiber distance, fiber shifting and fiber intercopenetration), it resulted that macroscopic porosity plays a key role dominating on all geometrical parameters; the global trend of stiffness and strength vs macroscopic porosity seems to follow one single third order polynomial law. Furthermore, for a given porosity, the fiber shifting between two layers affects the elasticity and strength substantially. In particular, the results shown in Fig. 5 indicate that approximately 30% of strength reduction is achieved when shifting the fiber to its maximum value (α = 0.5) .
More fundamental studies on cellular solids [27] have shown that the elastic properties of an open cell structure exhibits a power law relationship with porosity in the form E = E s (1 − φ) n in which E s is the intrinsic elastic modulus, while the exponent n is in the order of 1 to 3 depending on the deformation mechanisms occurring upon loading. For structures mostly subjected to axial loading n tends to 1; while for the bending dominated microstructures the exponent tends to 2. Experimental data on trabecular bone reported in Ref. [28] also show macroscopic elastic modulus exhibiting a power law dependence with n = 3.
In our study, compressive strength is also ruled by a power law with n=3. Experimental data on trabecular bone reported in Ref. [28] also show macroscopic compressive strength which shows a power law dependence with n=2. This is motivated by the buckling dominated failure mode of thin trabeculae subjected to bending and axial load. The defor- mation mechanisms occurring in the 3D printed scaffolds under study are mainly compression on the fiber intersection for scaffolds with α = 0, shear stress for small values of α and mainly bending without axial load for α = 1/2 (i.e. the maximum allowed fiber shifting). This deformation mechanisms may motivate the power law found in this study. In Ref. [13] glass scaffolds obtained through a different manufacturing technique are studied. The compressive and tensile strength of the CEL2 material obtained through the sponge replication process is analysed through laboratory tests and theoretical models. The microstructure of this particular scaffold does not exhibit a regular geometry like 3D printed scaffolds, therefore, a clear deformation mechanisms cannot be identified. For the sake of comparison with the results found in our study, the strength properties found in Ref. [13] for the CEL2 material under compression has been normalized by the theoretical strength of the material (C=33 MPa, see Table 3 in Ref. [13]). Figure 12 shows the normalized strength found in our study (filled bullet) and the normalized strength found in Ref. [13] (filled pentagons). Only data for porosity ranging approximatively φ = 0.55 − 0.75 are available; in this porosity   range the agreement with the results found with our scaffold is qualitatively satisfactory. It is worth to underline that, the comparison with the data from Ref. [13] are not meant to compare strength between the two different scaffolds, as the material and the microstructure are different. However, this agreement allows inferring that, in addition to the intrinsic strength of the constituent material, porosity is the major determinant in the compressive strength of glass scaffold. It remains to be determined whether this conclusion holds also for other loading conditions such as tensile or shear loading. A validation of the computational procedure has been carried out by comparing the computed strength with experimental results found in Ref. [6]. In this paper, robocast Sr-HT-Gahnite scaffolds have been characterized under compression. The geometrical features of the the scaffolds studied in Ref. [6] are the same as those considered in our study, although the constituent material is not the same. Under the assumption that the constituent material Sr-HT-Gahnite exhibits the same brittle nature as that of the glass scaffold object of this study, the comparison between the results of this study and that reported in Ref. [6] can be performed by selecting a suitable value of the intrinsic tensile strength σ t in the numerical simulation. The selection of the σ t parameter has been done by matching the macroscopic strength at a porosity of approximatively φ = 0.5 (a perfect match of the porosity was not achieved). The relationship between the macroscopic strength and the porosity, for porosity higher than φ = 0.5 was satisfactory (see Fig. 8). The proposed computational strategy features the typical limitations of simulations carried out using linear kinematic relationships (small strains and small displacements) in damaging elastic solids. Element annihilation as that used with voxel-like finite element models [29] implies that stress and strain fields in the annihilated elements are not reliable as large strains may locally occur as material undergoes the damage process. As the amount of incremental fracture surface is controlled (through a predetermined number of annihilated elements) and, as displacement is re-scaled at each step so to reduce stress within the strength limit, the stress and strain fields in the vicinity of the damaged elements can still be considered to be as accurate as those occurring in the original μ − CT based finite element models. The present computational approach is an incremental procedure that simulates neither a force-controlled, nor a displacement-controlled experiment. From the experimental point of view, it simulates a loading process with a fracture surface control. The inset in Fig. 2 shows a qualitative force displacement response of an elasto-damaging material; depending on the geometry and boundary conditions, this response may exhibit a snapback instability [30] (as that depicted in the inset), this phenomenon can be simulated through a fracture surface control simulation as the one reported in this paper. The expected experimental response carried out in displacement control would be the one repre-sented by the solid thick line reported in the inset of Fig. 2 obtained by considering a monotonically increasing displacement, this explains the oscillations occurring at the end of the stress vs strain plot in Fig. 10. The overall strength of the scaffold will be determined by two main parameters: i) the strength defined at the element length scale (σ t , σ c ), and ii) the number of elements to be removed at each increment N toll . Relationship between element size and intrinsic strength parameters should be also included on the basis of the typical flaw size occurring in a single element volume. XFEM simulations reported in Ref. [6] are obtained by using an intrinsic tensile strength of the constituent material of 97 MPa; this value is about 5 times lower than that used in the present work in the validation phase (Fig. 8); it is believed that this mismatch between the two intrinsic material properties is mainly due to the substantially smaller element size used in the present work with respect to the element size used in Refs. [6,31], although this latter can only be inferred through inspection of their Fig. 5 as it has not been explicitly reported in the paper. For the above reason, it must be emphasized that only an indirect validation was here provided as normalized values were compared and partial validation was achieved as few datapoints were available.
In general, higher strength for finer mesh should be allowed for as higher stress gradient are achieved in the proximity of defects. In principle, the above parameters should be identified through comparison with experimental tests performed on scaffold and relevant finite element modeling. In this study, the role of these parameters was not investigated as normalized results have been presented, limiting the validity of the study at a comparison between scaffold having different morphological parameters. The effect of the parameter N toll was studied by running three analyses in which three values of the parameter were used: N toll = 10 −5 N e , 1 2 × 10 −5 N e , 10 × 10 −5 N e . The results are shown in Fig. 13 in terms of normalized stress vs strain. The stresses are normalized with respect to the maximum stress obtained for the reference value N toll = 10 −5 N e . Figure 13 shows that increasing the number of elements to be removed at each step by a factor of 10 with respect to the reference value, the strength reduces of 4%; decreasing N toll to half of the reference value, the strength is increased of 9%, this should be taken into account when selecting the analyses parameters.
Furthermore, the role played by the real geometrical features as assessed through μ−CT scanning has been highlighted. Another major limitation on the presented results is owed to the limited size of the volume included in the analysis domain. Here 8 fiber intersections are included; models including larger volumes (with more fiber intersections) would have increased the chance to include localized defects of the printed structure. Therefore, the obtained normalized strength are to be considered as an upper limit of the scaffold strength, which is still an information of practical use.
Further studies should, therefore include: different loading mode (i.e. shear load) and/or loading directions, analysis of strength in different area of a scaffold (or to analyse larger volume) and eventually to determine the intrinsic material strength so to have the macroscopic strength instead of a relative (normalized) strength.

Conclusions
This work has presented a computational modeling for the prediction of the mechanical response of 3D printed bioactive glass scaffolds for bone tissue engineering. The finite element models are based on real three dimensional architecture of the printed scaffolds as obtained by μ − CT scanning of the samples. The model is expected to quantify the effect of the major morphological parameters on the elastic and strength properties of the scaffold allowing for possible imperfections resulting from the printing process. Selected design parameters have been investigated and the relevant impact on elasticity and strength has been highlighted. Through partial validation of the model we have shown that the proposed computational procedure is able to predict how mechanical properties of the scaffold change by changing geometrical parameters. This computational model provides a significant piece of information of practical use in the design process for bone tissue engineering scaffolds; this warrant further research work needed to overcome the limitations mentioned in the discussion section.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecomm ons.org/licenses/by/4.0/.