3D printing of porous scaffolds with controlled porosity and pore size values

materials Article 3D Printing of Porous Scaffolds with Controlled Porosity and Pore Size Values Irene Buj-Corral * ID , Ali Bagheri and Oriol Petit-Rojo Department of Mechanical Engineering, Universitat Politècnica de Catalunya, Av. Diagonal, 647, 08028 Barcelona, Spain; [email protected] (A.B.); [email protected] (O.P.-R.) *Correspondence: [email protected]; Tel.: +34-93-4054015 Received: 27 June 2018; Accepted: 22 August 2018; Published: 25 August 2018   Abstract: 3D printed scaffolds can be used, for example, in medical applications for simulating body tissues or for manufacturing prostheses. However, it is difficult to print porous structures of specific porosity and pore size values with fused deposition modelling (FDM) technology. The present paper provides a methodology to design porous structures to be printed. First, a model is defined with some theoretical parallel planes, which are bounded within a geometrical figure, for example a disk. Each plane has randomly distributed points on it. Then, the points are joined with lines. Finally, the lines are given a certain volume and the structure is obtained. The porosity of the structure depends on three geometrical variables: the distance between parallel layers, the number of columns on each layer and the radius of the columns. In order to obtain mathematical models to relate the variables with three responses, the porosity, the mean of pore diameter and the variance of pore diameter of the structures, design of experiments with three-level factorial analysis was used. Finally, multiobjective optimization was carried out by means of the desirability function method. In order to favour fixation of the structures by osseointegration, porosity range between 0.5 and 0.75, mean of pore size between 0.1 and 0.3 mm, and variance of pore size between 0.000 and 0.010 mm 2 were selected. Results showed that the optimal solution consists of a structure with a height between layers of 0.72 mm, 3.65 points per mm 2 and a radius of 0.15 mm. It was observed that, given fixed height and radius values, the three responses decrease with the number of points per surface unit. The increase of the radius of the columns implies the decrease of the porosity and of the mean of pore size. The decrease of the height between layers leads to a sharper decrease of both the porosity and the mean of pore size. In order to compare calculated and experimental values, scaffolds were printed in polylactic acid (PLA) with FDM technology. Porosity and pore size were measured with X-ray tomography. Average value of measured porosity was 0.594, while calculated porosity was 0.537. Average value of measured mean of pore size was 0.372 mm, while calculated value was 0.434 mm. Average value of variance of pore size was 0.048 mm 2 , higher than the calculated one of 0.008 mm 2 . In addition, both round and elongated pores were observed in the printed structures. The current methodology allows designing structures with different requirements for porosity and pore size. In addition, it can be applied to other responses. It will be very useful in medical applications such as the simulation of body tissues or the manufacture of prostheses. Keywords: fused deposition modeling; 3D printing; scaffolds; porosity; pore size; multiobjective optimization 1. Introduction Many new applications have arisen as a result of recent advances in 3D printing techniques. For example, printed parts are used for manufacturing space instrumentation, for both prototypes and flying parts [ 1 ], for manufacturing cost-effective parts in the sports industry or for developing new Materials 2018,11, 1532; doi:10.3390/ma11091532 www.mdpi.com/journal/materials Materials 2018,11, 1532 2 of 18 protective structures for vehicles in the automotive industry [ 2 ]. 3D printing has many different medical applications, such as bioprinting tissues and organs, building vascularized organs, the manufacture of customized implants, prostheses and models for surgical preparation, among others [ 3 ]. Thus, 3D printed scaffolds can be employed as templates for initial cell attachment and tissue formation, for example in bone tissue engineering [ 4 ]. They can also be used for fixing prostheses by means of osseointegration. Scaffolds could be printed in different materials such as titanium, degradable polymers and degradable ceramics [ 5 , 6 ]. Specifically, in FDM or fused deposition modelling technique a material is melted through an extrusion head and deposited layer by layer [ 7 ]. Main advantages of FDM technology are its easiness of use and the fact that it allows printing a wide range of materials, as long as they can be extruded, for example plastic materials such as acrylonitrile butadiene styrene (ABS) or polylactic acid (PLA). It is also more cost-effective than other additive manufacturing techniques, and its lead times are short. However, it also has disadvantages. It does not provide high dimensional precision, layer steps are usually observed on the part’s surface, causing the surfaces not to be smooth. Besides, the use of and the use of scaffolds with biocompatible materials is difficult, because the technology is limited to materials whose viscosity is sufficiently low that they can be extruded, but high enough so that their shape is maintained after extrusion [8–10]. Regarding the design of 3D printed porous scaffolds that simulate tissues, some properties to keep in mind are: surface area and interconnectivity, which are related to cell growth; permeability, which governs nutrient transport; and mechanical strength, which assures support and protection, among other properties. One possibility to achieve required porosity is to use hierarchical scaffold design, creating libraries of unit cells that can be joined to obtain scaffold structures. Hollister observed that increasing the material volume of a certain structure increased elastic modulus and decreased permeability [ 11 ]. In the same line, Egan et al. defined four different types of scaffold structures, taking into account either beam-based unit cells or truss-based unit cells. In addition, each structure could be created by means of either continuous or hierarchical patterning. They found that, given a certain porosity value, truss-based scaffolds have higher surface areas and lower elastic moduli than beam-based ones [ 12 ]. On the other hand, Arabnejad et al. presented a visualization method that allows understanding the relationship among cell topology, pore size and porosity in stretch-dominated structures, such as tetrahedron and octet trusses [13]. Another possibility to create structures with required porosity is use of topology optimization, which consists of distributing material in regions having low and high material density respectively. This is achieved by periodically repeating a unit cell, which is composed of areas with and without material [ 14 , 15 ]. When applying topological optimization methods, several authors have carried out multiobjective optimization of different properties in scaffolds. Lin et al. used an objective function that assigned different weights to the two responses considered: porosity and stiffness. The function was then maximized or minimized by numerical methods [ 16 ]. In a similar way, Guest and Prévost [ 17 ], as well as Hollister and Lin [ 18 ] used topology optimization to maximize stiffness and fluid permeability. Kang et al. used a homogenization-based topology optimization method to achieve the required bulk modulus and isotropic diffusivity for a certain porosity value [ 19 ]. Other properties of the structures, such as thermal conductivity [ 20 ], have been addressed. On the other hand, although the increase in surface area of pores helps tissue growth, this growth is facilitated if concave surfaces are used. In this direction, Egan et al. modelled scaffolds with curvature. They first fixed required porosity to assure a certain permeability value, which is desirable for nutrient transport, and then addressed the problem of mechanical strength of structures [21]. Different authors have tested the properties of FDM printed porous scaffolds. Regarding mechanical properties, Habib et al. used finite element analysis and compressive strength tests [ 22 ]. Wang et al. printed scaffolds for vascularized bone tissues in different materials, in order to test both their biomimetics and their strength [ 23 ]. Aw et al. tested the effect of printing parameters on tensile strength of conductive acrylonitrile butadiene styrene/zinc oxide (CABS/ZnO) composites [ 24 ]. Materials 2018,11, 1532 3 of 18 Helguero et al. modelled artificial bones and printed them in acrylonitrile butadiene styrene (ABS). They tested both anisotropy and compressive strength of the scaffolds [ 25 ]. As for porosity, Gregor et al. printed PLA scaffolds and measured their porosity by means of X-ray microtomography [ 26 ]. Regarding surface finish, Townsend et al. listed most usual methods for measuring roughness profiles (contact stylus), and surface topography (confocal microscopy, focus variation microscopy, coherence scanning interferometry, chromatic confocal microscopy, conoscopic holography, atomic force microscopy (AFM), and elastomeric sensors [ 27 ]). Krolczyk et al. compared the roughness obtained in turning processes with that obtained in FDM processes [ 28 ]. They observed that the machined surface had an anisotropic and periodic structure, while the printed surface had an undirected structure. With the manufacturing conditions employed, the FDM process showed higher roughness values than the turning process. In the present paper, a model was developed to define the pore size and the porosity of porous structures. Unlike other methods that are based on truss structures, the present model allows obtaining irregular porous structures from random location of columns in the space, which leave voids among them. Specifically, the structure was modelled with parallel planes joined by columns, with a certain number of columns on each plane. The model was applied to a disk shape. Three variables were defined: the distance between parallel planes, the number of base points for columns on each plane, and the radius of each column. Next, dimensional analysis was used to reduce the number of process variables to 2. Then, the requirements were defined for a specific application case: the use of a porous structure in external layers of hemispherical hip prostheses. Subsequently, the design of experiments, with three-level factorial analysis, was used to obtain mathematical models for porosity, mean of pore size and variance of pore size as a function of dimensionless variables. They allowed multiobjective optimization in order to determine the optimal values for the process parameters. In order to compare experimental results to computationally calculated ones, samples were printed with FDM technology, and total porosity, as well as pore size, was measured. X-ray tomography was used to determine the total porosity of the printed structures by means of computation of plastic volume and comparison with total volume of the printed shape. The present study will help designing and manufacturing porous structures with specific requirements regarding the porosity and pore size that favour osseointegration. The same methodology can be used, however, to achieve other requirements of porous structures, such as mechanical strength and/or nutrient transport. 2. Materials and Methods 2.1. Model for the Porous Structure 2.1.1. Model Definition The printing process of a random porous structure, such as a trabecular one, presents some difficulties: the walls in certain parts of the structure are too thin (Figure 1), there are some areas with burrs, some parts of the structure have a high inclination angle that leads to the use of printing structure supports, etc. Such difficulties can be attributed to the fact that, when designing the structure, a completely random distribution of points in space is used. Since distribution is random, connections between different points in space can have any orientation. In the present paper, a proposed solution for modelling porous structures is presented. It involves defining some theoretical parallel planes, each one of them with several theoretical points that are randomly distributed. Once the points have been created on the surfaces, it is necessary to connect them by means of theoretical lines. For doing this, each point is connected to the three nearest points of the same plane or of the plane that is immediately below. Thus, almost vertical lines will be created that will have a correct inclination angle for printing. Figure 2shows a scheme of the model applied to a disk shape. Materials 2018,11, 1532 4 of 18 Materials 2018, 11, x FOR PEER REVIEW 3 of 19 et al. printed PLA scaffolds and measured their porosity by means of X-ray microtomography [26]. Regarding surface finish, Townsend et al. listed most usual methods for measuring roughness profiles (contact stylus), and surface topography (confocal microscopy, focus variation microscopy, coherence scanning interferometry, chromatic confocal microscopy, conoscopic holography, atomic force microscopy (AFM), and elastomeric sensors [27]). Krolczyk et al. compared the roughness obtained in turning processes with that obtained in FDM processes [28]. They observed that the machined surface had an anisotropic and periodic structure, while the printed surface had an undirected structure. With the manufacturing conditions employed, the FDM process showed higher roughness values than the turning process. In the present paper, a model was developed to define the pore size and the porosity of porous structures. Unlike other methods that are based on truss structures, the present model allows obtaining irregular porous structures from random location of columns in the space, which leave voids among them. Specifically, the structure was modelled with parallel planes joined by columns, with a certain number of columns on each plane. The model was applied to a disk shape. Three variables were defined: the distance between parallel planes, the number of base points for columns on each plane, and the radius of each column. Next, dimensional analysis was used to reduce the number of process variables to 2. Then, the requirements were defined for a specific application case: the use of a porous structure in external layers of hemispherical hip prostheses. Subsequently, the design of experiments, with three-level factorial analysis, was used to obtain mathematical models for porosity, mean of pore size and variance of pore size as a function of dimensionless variables. They allowed multiobjective optimization in order to determine the optimal values for the process parameters. In order to compare experimental results to computationally calculated ones, samples were printed with FDM technology, and total porosity, as well as pore size, was measured. X-ray tomography was used to determine the total porosity of the printed structures by means of computation of plastic volume and comparison with total volume of the printed shape. The present study will help designing and manufacturing porous structures with specific requirements regarding the porosity and pore size that favour osseointegration. The same methodology can be used, however, to achieve other requirements of porous structures, such as mechanical strength and/or nutrient transport. 2. Materials and Methods 2.1. Model for the Porous Structure 2.1.1. Model Definition The printing process of a random porous structure, such as a trabecular one, presents some difficulties: the walls in certain parts of the structure are too thin (Figure 1), there are some areas with burrs, some parts of the structure have a high inclination angle that leads to the use of printing structure supports, etc. Figure 1. Trabecular structure with thin walls. Figure 1. Trabecular structure with thin walls. Materials 2018, 11, x FOR PEER REVIEW 4 of 19 Such difficulties can be attributed to the fact that, when designing the structure, a completely random distribution of points in space is used. Since distribution is random, connections between different points in space can have any orientation. In the present paper, a proposed solution for modelling porous structures is presented. It involves defining some theoretical parallel planes, each one of them with several theoretical points that are randomly distributed. Once the points have been created on the surfaces, it is necessary to connect them by means of theoretical lines. For doing this, each point is connected to the three nearest points of the same plane or of the plane that is immediately below. Thus, almost vertical lines will be created that will have a correct inclination angle for printing. Figure 2 shows a scheme of the model applied to a disk shape. Figure 2. Structure obtained after applying the model to a disk shape. The marching cubes algorithm allows lines and points to be given volume [29]. It is used in the present paper to create columns of the structure from the theoretical lines. The structure obtained after applying the model and the marching cubes algorithm to a disk shape is shown in Figure 3. Figure 3. Structure obtained after applying the marching cubes algorithm to a disk shape. Thus, there are three parameters to be controlled: distance between parallel planes H (mm), number of points per surface unit on a certain plane N (points/mm2) (Figure 4a), and radius r of columns that are created around the model lines (mm) (Figure 4b). Figure 4. Parameters that define the porous structures: (a) parameters H and N, and (b) parameter r. With the aim of determining optimum values for parameters, it is necessary to find the relationship between the three model parameters and the responses of porosity and pore size. To do so, a design of experiments was used to obtain mathematical models and perform multiobjective a b Figure 2. Structure obtained after applying the model to a disk shape. The marching cubes algorithm allows lines and points to be given volume [ 29 ]. It is used in the present paper to create columns of the structure from the theoretical lines. The structure obtained after applying the model and the marching cubes algorithm to a disk shape is shown in Figure 3. Materials 2018, 11, x FOR PEER REVIEW 4 of 19 Such difficulties can be attributed to the fact that, when designing the structure, a completely random distribution of points in space is used. Since distribution is random, connections between different points in space can have any orientation. In the present paper, a proposed solution for modelling porous structures is presented. It involves defining some theoretical parallel planes, each one of them with several theoretical points that are randomly distributed. Once the points have been created on the surfaces, it is necessary to connect them by means of theoretical lines. For doing this, each point is connected to the three nearest points of the same plane or of the plane that is immediately below. Thus, almost vertical lines will be created that will have a correct inclination angle for printing. Figure 2 shows a scheme of the model applied to a disk shape. Figure 2. Structure obtained after applying the model to a disk shape. The marching cubes algorithm allows lines and points to be given volume [29]. It is used in the present paper to create columns of the structure from the theoretical lines. The structure obtained after applying the model and the marching cubes algorithm to a disk shape is shown in Figure 3. Figure 3. Structure obtained after applying the marching cubes algorithm to a disk shape. Thus, there are three parameters to be controlled: distance between parallel planes H (mm), number of points per surface unit on a certain plane N (points/mm2) (Figure 4a), and radius r of columns that are created around the model lines (mm) (Figure 4b). Figure 4. Parameters that define the porous structures: (a) parameters H and N, and (b) parameter r. With the aim of determining optimum values for parameters, it is necessary to find the relationship between the three model parameters and the responses of porosity and pore size. To do so, a design of experiments was used to obtain mathematical models and perform multiobjective a b Figure 3. Structure obtained after applying the marching cubes algorithm to a disk shape. Thus, there are three parameters to be controlled: distance between parallel planes H (mm), number of points per surface unit on a certain plane N (points/mm 2 ) (Figure 4a), and radius r of columns that are created around the model lines (mm) (Figure 4b). Materials 2018, 11, x FOR PEER REVIEW 4 of 19 Such difficulties can be attributed to the fact that, when designing the structure, a completely random distribution of points in space is used. Since distribution is random, connections between different points in space can have any orientation. In the present paper, a proposed solution for modelling porous structures is presented. It involves defining some theoretical parallel planes, each one of them with several theoretical points that are randomly distributed. Once the points have been created on the surfaces, it is necessary to connect them by means of theoretical lines. For doing this, each point is connected to the three nearest points of the same plane or of the plane that is immediately below. Thus, almost vertical lines will be created that will have a correct inclination angle for printing. Figure 2 shows a scheme of the model applied to a disk shape. Figure 2. Structure obtained after applying the model to a disk shape. The marching cubes algorithm allows lines and points to be given volume [29]. It is used in the present paper to create columns of the structure from the theoretical lines. The structure obtained after applying the model and the marching cubes algorithm to a disk shape is shown in Figure 3. Figure 3. Structure obtained after applying the marching cubes algorithm to a disk shape. Thus, there are three parameters to be controlled: distance between parallel planes H (mm), number of points per surface unit on a certain plane N (points/mm2) (Figure 4a), and radius r of columns that are created around the model lines (mm) (Figure 4b). Figure 4. Parameters that define the porous structures: (a) parameters H and N, and (b) parameter r. With the aim of determining optimum values for parameters, it is necessary to find the relationship between the three model parameters and the responses of porosity and pore size. To do so, a design of experiments was used to obtain mathematical models and perform multiobjective a b Figure 4. Parameters that define the porous structures: (a) parameters H and N, and (b) parameter r. Materials 2018,11, 1532 5 of 18 With the aim of determining optimum values for parameters, it is necessary to find the relationship between the three model parameters and the responses of porosity and pore size. To do so, a design of experiments was used to obtain mathematical models and perform multiobjective optimization. Before optimization, a reduction of the variables was applied so that each response considered depends on two dimensionless variables. 2.1.2. Reduction of Variables Dimensional analysis allows reducing number of variables in a certain system [ 30 ]. If the problem depends on n dimensional variables, dimensional analysis reduces the problem to k dimensionless variables, where reduction is n − k = 1, 2, 3, or 4 depending on the complexity of the system. One of the most used techniques for reducing the number of variables is the Buckingham Pi theorem. This allows obtaining dimensionless parameters that are a product of powers named Π1,Π2,Π3, etc. [31]. First part of the Pi theorem explains how to find reduction j Equation (1): j=n−k (1) where n is the number of variables of the problem and k is the number of dependent variables. From n variables, k variables are chosen as basic and the rest are scale variables. The second part of the Pi theorem shows that each dimensionless parameter is defined as the product of the rest of variables with a proper exponent different than zero. In this way, all dimensionless groups are independent. Each dimensionless parameter is related to the other ones by means of a function Equation (2): Π1=g(Π2,Π3,Π4, . . .)(2) where Π1is a dimensionless variable, and Π2,Π3,Π4are other variables. In the present paper, three responses were analysed separately: porosity p (dimensionless), mean of pore size M (length) and variance of pore size V (length 2 ), with process variables H (length), N (length-2) and r (length). Variable r was chosen as the scale variable. For porosity, dimensionless parameters are described in Equations (3)–(5): Π1= (H[long1])1·(r[long1])−1=H r[dimensionless](3) Π2= (N[long−2])1·(r[long1])2=N·r2[dimensionless](4) Π3=p[dimensionless](5) and the function to be found is presented in Equation (6): Π3=g(Π1,Π2)→p=gH r, N·r2(6) Similar expressions were found for average value (M) and variance (V) of the probability density function and of the pore diameters Equations (7) and (8): M r=jH r, N·r2(7) V r2=kH r, N·r2(8) Thus, it was only necessary to perform experiments with two dimensionless factors, H/r and N·r2. Materials 2018,11, 1532 6 of 18 2.1.3. Requirements of the Structure The kind of porous structure that is defined in the present paper can be used, for example, to replace or simulate bone trabecular tissue. At the macroscopic level, bones are formed by cortical and trabecular compartments. Bone tissues have several functions, for example, support, protection, mineral storage, and nutrient transport [ 32 ]. If a printed structure is to be fixed by means of osseointegration, a certain degree of porosity is required. Two main factors governing the porous external surface are the pore size and the porosity of the structure. According to Karageorgiou [ 33 ] and Baino [ 34 ], optimum values of pore size were considered to range between 100 and 500 µ m for a trabecular structure. The mean porosity value of scaffolds is recommended to be between 50 and 75% in volume. Although a certain variability of pore size, with small pores that improve cell attachment and large pores that favour nutrient transport could be desirable in some cases [ 35 ], in the present work only osseointegration was considered. For this reason, the variance of the pore size should take the lowest possible value, in order to assure that the maximum number of pores would lie within the interval. High connectivity between pores in the structures is also required, as porosity has to be accessible; an interconnecting channel structure is needed to allow tissue to grow on it. 2.2. Design of Experiments In order to determine the required values for the two model variables considered and to optimize porosity, design of experiments was performed with simulated experiments. A full factorial design was employed with two dimensionless factors: H/r, the distance between surfaces where points are generated over a radius of generated volumes, and N · r 2 , the number of points per area unit. Three levels were defined for each factor and three replicates were performed for each experiment (since variability among replicates is high). In all, 3 × 3 2 = 27 runs were carried out. The radius of columns remained constant in all experiments. Levels for variables were selected according to previous tests. A summary of levels is presented in Table 1. Table 1. Levels of the factorial design. Level 1 Level 2 Level 3 H/r 3.2 4.0 4.8 N·r20.054 0.072 0.090 Responses considered were porosity, mean of pore diameters (mm), and variance of pore diameters (mm2), in order to consider both the mean and dispersion of the probability function of porosity. In order to calculate theoretical porosity and distribution of pore diameters, a program created by Dupuy et al. was used [ 36 ]. It is a MATLAB (MathWorks, Natick, MA, USA) script that, from a binary 3D image, returns the function of cumulated porosity vs. pore diameter. In the present study, a distribution of points in space was created with the program Rhinoceros 3D with the Grasshopper plug-in (Robert McNeel & Associates, Seattle, WA, USA). 2.3. Multiobjective Optimization Once the mathematical models had been obtained, it was necessary to choose values for the parameters (N, H, r) that provided the required values for porosity, average pore size, and variance, taking into account the requirements of the responses. The desirability function method was employed for the multiobjective optimization of responses [ 37 ]. It is based on defining desirability functions di(xi) for each response. The range for the desirability function is [0, 1]. Value 0 is assigned to a situation that is undesirable and value 1 is assigned to a completely desirable situation. Materials 2018,11, 1532 7 of 18 The objective function to maximize is the geometric average of the group of desirability functions. Taking into account different weights kifor each objective Equation (9): D=∑n i=1ki sn ∏ i=1 (di(xi))ki(9) where n is number of responses to be optimized at the same time, and ki defines the importance of each desirability function inside the objective function. There are three kinds of desirability functions depending on the optimization goal: maximizing, minimizing, and target value desirability functions. As an example, in Equation (10) the maximization function is defined: di,max(xi)=         0 xi≤L xi−L U−LbL<xi<U 1 U ≤xi (10) where L is a lower limit value from which the function becomes undesirable (di = 0) and U is an upper limit value from which the function is considered to be maximized (di = 1). Exponent b is a parameter that models the shape of the central part of the desirability function. It can take any value within [0, ∞ ]. If b takes values [0, <1], the central part of the function will be concave, and a small increment of the response from the unacceptable point leads to a great increment of d i . If b takes a value of 1, the central part of the function will be a straight line with constant slope. If b takes values [>1, ∞ ], the central part of the curve will be convex, and, until a value that is very similar to the optimal one is obtained, the d i value will be almost zero. Table 2shows the selected values of parameters for optimization. Table 2. Selected values of parameters for optimization. Response Type of Desirability Function Lower Value L Target Value T Upper Value U b Porosity Of target value 0.500 0.625 0.750 1 Mean (mm) Of target value 0.100 0.300 0.500 1 Variance (mm 2 ) Minimization 0.000 * 0.010 1 * means that a certain parameter was not used. For porosity and mean of pore size, lower (L) and upper (U) values were selected, defining an interval for each response. The average value of the interval was selected as target value (T). For variance of pore size, lower value (L) was 0, since variance is to be minimized so that the maximum number of pores is contained within the required interval for the pore size. The upper value (U) was set at 0.01. It is supposed that pore size follows a normal distribution, in which 95.4% of the population is included within the µ± 2 ·σ interval. By equalling the range of pore diameters [0.100, 0.500] with the a.m. interval Equation (11), the maximal standard deviation σmax = 0.100 mm is found, according to Equation (12). Thus, the maximal variance Variance max is 0.010 mm 2 according to Equation (13): [0.300 −2·σmax, 0.300 +2·σmax]mm = [0.100, 0.500]mm (11) 0.300 −0.100 2=0.100 mm =σmax (12) Variancemax =σmax2=0.010 mm2(13) A MATLAB program was implemented in a discrete optimization of the objective function by means of sweeping. The program generates mi equispaced values within an interval for each variable (H, N, and r). Then, it checks if each point (combination of values for the different variables), meets Materials 2018,11, 1532 8 of 18 the conditions defined in Table 1. Such conditions correspond to the space that is considered in the regression. Table 3corresponds to the intervals for variables and the m-values used for each variable. The value for m was selected so that the distance between contiguous values is lower than or equal to 0.005. Table 3. Intervals for variables and m-values used for optimization. Interval m H (mm) [0.48, 2.00] 300 N (points/mm2)[0.30, 4.00] 740 r (mm) [0.15, 0.41] 200 Different solutions to be found will depend on the importance ( ki ) of each response. Specifically, four different cases were analyzed: - The first case involves giving the same importance to all responses. - The second case consists in giving higher importance to porosity response over the rest of the responses (kporosity =5). - The third case gives a higher importance to mean (kmean =5). - The fourth case consists in giving higher importance to variance (kvariance =5) 2.4. Experimental Tests Printing experiments were performed in order to determine the experimental porosity of the samples. The techniques of fused filament fabrication (FFF) or fused deposition modelling (FDM) were selected. Cura software was used for generating the g-code for printing. Disk-shaped samples of 6 mm diameter and 3 mm height were designed. Selected parameters for the structure were H = 0.72 mm, r = 0.15 mm, and N = 3.65 points/mm 2 (according to results of multiobjective optimization in Section 3.2). Three specimens were printed in a dual-extruder Sigma printer from BCN3D. Due to the difficulty to print the required columns of diameter 0.3 mm, the geometry of the disks was rescaled by a scaling factor of five before printing (30 mm diameter and 15 mm height). The printing speed of the head was set to 37 mm/s. The nozzle diameter was 0.2 mm, the layer height was 0.1 mm, and the infill was 40%. The shell thickness was 0.4 mm. The geometry of the samples was measured by means of X-ray tomography with Zeiss Metrotom 800 equipment. From the geometry, it was possible to calculate the percentage of air (porosity) and the percentage of plastic in the structure with VG Studio Max software. From the 3D geometry, a cross-section was obtained at half the height of each specimen, where the width of pores was measured. 3. Results In the present section results about the design of the experiments, multiobjective optimization, and experimental tests are presented. 3.1. Design of Experiments Table 4shows the results for the 27 different runs performed. The levels for the two factors considered are presented in codified units. Results for porosity, the mean of the pore size, and the variance of the pore size were simulated with a MATLAB program. Since the results of experiments come from simulation, the order of the runs is indifferent, and a standard order was used for performing the experiments. Materials 2018,11, 1532 9 of 18 Table 4. Results of the porosity, the mean of the pore size, and the variance of the pore size. Nr H/r N·r2Porosity Mean (mm) Variance (mm2) 1−1−1 0.574 0.524 0.009 2−1 0 0.491 0.439 0.009 3−1 1 0.409 0.359 0.006 4 0 −1 0.599 0.549 0.012 5 0 0 0.537 0.450 0.011 6 0 1 0.472 0.378 0.005 7 1 −1 0.645 0.578 0.015 8 1 0 0.563 0.466 0.007 9 1 1 0.500 0.414 0.008 10 −1−1 0.587 0.559 0.022 11 −1 0 0.473 0.414 0.008 12 −1 1 0.436 0.369 0.007 13 0 −1 0.635 0.547 0.008 14 0 0 0.552 0.482 0.012 15 0 1 0.484 0.401 0.008 16 1 −1 0.644 0.575 0.010 17 1 0 0.582 0.491 0.012 18 1 1 0.521 0.411 0.006 19 −1−1 0.549 0.525 0.016 20 −1 0 0.477 0.447 0.010 21 −1 1 0.430 0.397 0.008 22 0 −1 0.621 0.542 0.009 23 0 0 0.516 0.439 0.006 24 0 1 0.462 0.376 0.006 25 1 −1 0.641 0.575 0.015 26 1 0 0.571 0.476 0.008 27 1 1 0.502 0.401 0.009 Data were analysed with Minitab 17 (Minitab, State College, PA, USA). Multiple linear regressions were used to examine the mathematical models for the porosity, the mean of pore size, and the variance of the pore size as a function of the two dimensionless factors considered. For all models, hypotheses were formulated for the linearity, constant variance, normality, and independence of errors. All terms having p-values higher than 5% were not included in the model. Equation (14) corresponds to porosity in real units: Porosity =0.614 +0.051H rreal −3.948(N·r2)real (14) The R-adj coefficient is 95.89%. The function for the mean pore size in real units is shown in Equation (15): Mean =0.929 +0.024H rreal −11.336(N·r2)real +47.222(N·r2)2 real (15) Materials 2018,11, 1532 16 of 18 height of steps or controlling the flow in order to avoid filaments within the pores [ 40 ]. Other additive manufacturing technologies, such as stereolithography (STL), having higher dimensional accuracy than FDM, can also be employed [41]. 5. Conclusions In the present paper a new method for designing porous structures is presented. First, a geometric model is defined that is based on parallel layers with several columns of a certain radius that join them. The model can be applied to a certain shape, for example a disk. Then, by means of dimensional analysis, number of variables is reduced. In this case, two dimensionless variables are used: H/r and N · r 2 . Next, design of experiments allows obtaining mathematical models that relate process variables to responses, in this case the porosity (dimensionless), the mean of the pore size (mm) and the variance of the pore size (mm 2 ) of the structures. Then, multiobjective optimization provides the values of the variables that allow obtaining required values for the responses. If the model is applied to a trabecular structure, in order to obtain target porosity of 0.625, target mean of pore size of 0.300 mm and low variance of pore size, recommended values were found for model variables: lowest radius of columns r = 0.15 mm and a medium value for height between parallel surfaces of H = 0.72. The number of points ranges from N = 3.20 to 3.90 points/mm 2 depending on the importance given to each response. For example, if the three responses have the same importance, N = 3.65 points/mm 2 are recommended. Higher number of points is recommended when mean or variance of pore size is more important than the rest of the variables. On the contrary, lower number of points is recommended when porosity is more important than the rest of the variables. The porosity and the mean of pore size obtained in experimental tests printed in polylactic acid (PLA), 0.594 and 0.372 mm respectively, were in good agreement with the simulated porosity and mean of pore size of 0.537 and 0.434 mm, respectively. Experimental variance of pore size is higher than the calculated one. For fixed values of the distance between layers and of the radius of columns, the greater the number of points on layers, the lower the calculated porosity and the mean of pore size. A larger number of columns means there is less room for pores, leading to lower porosity, as well as shorter distance among columns, with lower pore size. The variance of pore size also decreases with number of points per surface unit. By varying the distance between layers, the radius of the columns and the number of columns per unit area, structures with different porosity and pore size can be achieved. The higher distance between layers is the higher porosity and mean of pore size. The higher the radius of columns is the lower porosity and mean of pore size. Thus, the model will be very useful in defining the required values for geometric variables, which will allow a certain degree of porosity and pore size in porous structures to be obtained. The main limitation of the methodology is that it only takes into account the porosity of the structures that is required to achieve osseointegration, but not their mechanical properties (elastic modulus, compressive strength) or their ability to transport nutrients (permeability, diffusivity). Another limitation is that the structures are difficult to print with FDM technologies because of the low diameter of columns which is only slightly higher than nozzle diameter employed. For this reason, in the present work geometry was rescaled by a scaling factor of five prior to printing. Surface finish presents the stair-stepping defect due to layer-by-layer deposition. In future work, other requirements for structures, related to either mechanical strength or mass transport, will be addressed. In addition, improvement of the FDM printing process is required in order to obtain more accurate and smooth parts. For example, nozzle diameter can be reduced and nozzle retraction can be controlled in order to minimise the stair-stepping effect and the presence of threads within pores respectively. Author Contributions: I.B.-C. conceptualized the research, wrote the original draft and reviewd it, A.B. performed data analysis and calculated results, O.P.-R. defined the methodology. Funding: This research was funded by Spanish Ministry of Economy, Industry, and Competitiveness grant number DPI2016-80345R. Materials 2018,11, 1532 17 of 18 Acknowledgments: The authors thank Ramón Casado-López and Sergi Aymerich for their support with experimental tests. Conflicts of Interest: The authors declare no conflicts of interest. References 1. Rochus, P.; Plesseria, J.-Y.; Van Elsen, M.; Kruth, J.-P.; Carrus, R.; Dormal, T. New applications of rapid prototyping and rapid manufacturing (RP/RM) technologies for space instrumentation. Acta Astronaut. 2007,61, 352–359. [CrossRef] 2. Ngo, T.D.; Kashani, A.; Imbalzano, G.; Nguyen, K.T.Q.; Hui, D. Additive manufacturing (3D printing): A review of materials, methods, applications and challenges. Compos. Part B Eng. 2018 ,143, 172–196. [CrossRef] 3. Ventola, C.L. Medical Applications for 3D Printing: Current and Projected Uses. Pharm. Ther. 2014 ,39, 704–711. 4. Leukers, B.; Gülkan, H.; Irsen, S.H.; Milz, S.; Tille, C.; Schieker, M.; Seitz, H. Hydroxyapatite scaffolds for bone tissue engineering made by 3D printing. J. Mater. Sci. Mater. Med. 2005 ,16, 1121–1124. [CrossRef] [PubMed] 5. Hollister, S.J.; Lin, C.Y.; Saito, E.; Lin, C.Y.; Schek, R.D.; Taboas, J.M.; Williams, J.M.; Partee, B.; Flanagan, C.L.; Diggs, A.; et al. Engineering craniofacial scaffolds. Orthod. Craniofac. Res. 2005 ,8, 162–173. [CrossRef] [PubMed] 6. Klammert, U.; Vorndran, E.; Reuther, T.; Müller, F.A.; Zorn, K.; Gbureck, U. Low temperature fabrication of magnesium phosphate cement scaffolds by 3D powder printing. J. Mater. Sci. Mater. Med. 2010 ,21, 2947–2953. [CrossRef] [PubMed] 7. An, J.; Teoh, J.E.M.; Suntornnond, R.; Chua, C.K. Design and 3D Printing of Scaffolds and Tissues. Engineering 2015,1, 261–268. [CrossRef] 8. Lipson, H.; Kurman, M. Fabricated: The New World of 3D Printing; Wiley and Sons: New York, NY, USA, 2013. 9. Gibson, I.; Rosen, D.; Stucker, B. Additive Manufacturing Technologie. Springer: New York, NY, USA, 2015. 10. Chia, H.N.; Wu, B.M. Recent advances in 3D printing of biomaterials. J. Biol. Eng. 2015 ,9, 4. [CrossRef] [PubMed] 11. Hollister, S.J. Porous scaffold design for tissue engineering. Nat. Mater. 2005 ,4, 518–524. [CrossRef] [PubMed] 12. Egan, P.F.; Ferguson, S.J.; Shea, K. Design of hierarchical three-dimensional printed scaffolds considering mechanical and biological factors for bone tissue engineering. J. Mech. Des. N. Y. 2017 ,139, 061401-7. [CrossRef] 13. Arabnejad, S.; Burnett Johnston, R.; Pura, J.A.; Singh, B.; Tanzer, M.; Pasini, D. High-strength porous biomaterials for bone replacement: A strategy to assess the interplay between cell morphology, mechanical properties, bone ingrowth and manufacturing constraints. Acta Biomater. 2016 ,30, 345–356. [CrossRef] [PubMed] 14. Bendsøe, M.P.; Kikuchi, N. Generating optimal topologies in structural design using a homogenization method. Comput. Methods Appl. Mech. Eng. 1988,71, 197–224. [CrossRef] 15. Sigmund, O. Materials with prescribed constitutive parameters: An inverse homogenization problem. Int. J. Solids Struct. 1994,31, 2313–2329. [CrossRef] 16. Lin, C.Y.; Kikuchi, N.; Hollister, S.J. A novel method for biomaterial scaffold internal architecture design to match bone elastic properties with desired porosity. J. Biomech. 2004,37, 623–636. [CrossRef] [PubMed] 17. Guest, J.K.; Prévost, J.H. Optimizing multifunctional materials: Design of microstructures for maximized stiffness and fluid permeability. Int. J. Solids Struct. 2006,43, 7028–7047. [CrossRef] 18. Hollister, S.J.; Lin, C.Y. Computational design of tissue engineering scaffolds. Comput. Methods Appl. Mech. Eng. 2007,196, 2991–2998. [CrossRef] 19. Kang, H.; Lin, C.-Y.; Hollister, S.J. Topology optimization of three dimensional tissue engineering scaffold architectures for prescribed bulk modulus and diffusivity. Struct. Multidiscip. Optim. 2010 ,42, 633–644. [CrossRef] 20. De Kruijf, N.; Zhou, S.; Li, Q.; Mai, Y.W. Topological design of structures and composite materials with multiobjectives. Int. J. Solids Struct. 2007,44, 7092–7109. [CrossRef] Materials 2018,11, 1532 18 of 18 21. Egan, P.F.; Shea, K.A.; Ferguson, S.J. Simulated tissue growth for 3D printed scaffolds. Biomech. Model. Mechanobiol. 2018. [CrossRef] [PubMed] 22. Habib, F.N.; Nikzad, M.; Masood, S.H.; Saifullah, A.B.M. Design and Development of Scaffolds for Tissue Engineering Using Three-Dimensional Printing for Bio-Based Applications. 3D Print. Addit. Manuf. 2016 ,3, 119–127. [CrossRef] 23. Wang, M.O.; Vorwald, C.E.; Dreher, M.L.; Mott, E.J.; Cheng, M.-H.; Cinar, A.; Mehdizadeh, H.; Somo, S.; Dean, D.; Brey, E.M.; et al. Evaluating 3D-Printed Biomaterials as Scaffolds for Vascularized Bone Tissue Engineering. Adv Mater. 2015,27, 138–144. [CrossRef] [PubMed] 24. Aw, Y.Y.; Yeoh, Ch.K.; Idris, M.A.; Teh, P.L.; Hamzah, K.A. Effect of printing parameters on tensile, dynamic mechanical, and thermoelectric properties of FDM 3D printed CABS/ZnO composites. Materials 2018 , 11, 466. [CrossRef] [PubMed] 25. Helguero, C.G.; Amaya, J.L.; Komatsu, D.E.; Pentyala, S.; Mustahsan, V.; Ramirez, E.A.; Kao, I. Trabecular Scaffolds’ Mechanical Properties of Bone Reconstruction Using Biomimetic Implants. Procedia CIRP 2017 ,65, 121–126. [CrossRef] 26. Gregor, A.; Filová, E.; Novák, M.; Kronek, J.; Chlup, H.; Buzgo, M.; Blahnová, V.; Lukášová, V.; Bartoš, M.; Neˇcas, A.; et al. Designing of PLA scaffolds for bone tissue replacement fabricated by ordinary commercial 3D printer. J. Biol. Eng. 2017,11, 31. [CrossRef] [PubMed] 27. Townsend, A.; Senin, N.; Blunt, L.; Leach, R.K.; Taylor, J.S. Surface texture metrology for metal additive manufacturing: A review. Precis. Eng. 2016,46, 34–47. [CrossRef] 28. Krolczyk, G.; Raos, P.; Legutko, S. Experimental Analysis of Surface Roughness and Surface Texture of Machined and Fused Deposition Modelled Parts. Tehniˇcki Vjesnik 2014,21, 217–221. 29. Lorensen, W.E.; Cline, H.E. Marching cubes: A high resolution 3D surface construction algorithm. In Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques—SIGGRAPH ’87; ACM Press: New York, NY, USA, 1987; pp. 163–169. [CrossRef] 30. White, F.M. Fluid Mechanics; McGraw Hill: New York, NY, USA, 2011. 31. Buckingham, E. On Physically Similar Systems: Illustrations of the Use of Dimensional Equations. Phys. Rev. 1914,4, 345–376. [CrossRef] 32. De Acevedo-Marques, P.M.; Mencattini, A.; Salmeri, M.; Rangayyan, R.M. Medical Image Analysis and Informatics: Computer-Aided Diagnosis and Therapy; CRC Press: Boca Raton, FL, USA, 2018. 33. Karageorgiou, V.; Kaplan, D. Porosity of 3D biomaterial scaffolds and osteogenesis. Biomaterials 2005 ,26, 5474–5491. [CrossRef] [PubMed] 34. Baino, F.; Vitale-Brovarone, C. Mechanical properties and reliability of glass-ceramic foam scaffolds for, bone repair. Mater. Lett. 2014,118, 27–30. [CrossRef] 35. Egan, P.F.; Gonella, V.C.; Engensperger, M.; Ferguson, S.J.; Shea, K. Computationally designed lattices with tuned properties for tissue engineering using 3D printing. PLoS ONE 2017,12, 1–21. [CrossRef] [PubMed] 36. Dupuy, P.M.; Austin, P.; Delaney, G.W.; Schwarz, M.P. Pore scale definition and computation from tomography data. Comput. Phys. Commun. 2011,182, 2249–2258. [CrossRef] 37. Derringer, G.C.; Suich, R. Simultaneous optimization of several response variables. J. Qual. Technol. 1980 ,12, 214–219. [CrossRef] 38. Rangaiah, G.P. Multi-Objective Optimization: Techniques and Applications in Chemical Engineering; World Scientific: Singapore, 2009. 39. Zhao, H.; Li, L.; Ding, S.; Liu, C.; Ai, J. Effect of porous structure and pore size on mechanical strength of 3D-printed comby scaffolds. Mater. Lett. 2018,223, 21–24. [CrossRef] 40. Popescu, D.; Zapciu, A.; Amza, C.; Baciu, F.; Marinescu, R. FDM process parameters influence over the mechanical properties of polymer specimens: A review. Polym. Test. 2018,69, 157–166. [CrossRef] 41. Melchels, F.P.W.; Feijen, J.; Grijpma, D.W. A review on stereolithography and its applicationsin biomedical engineering. Biomaterials 2010,31, 6121–6130. [CrossRef] [PubMed] © 2018 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).