Prediction of permeability of regular scaffolds for skeletal tissue engineering: A combined computational and experimental study

Acta Biomaterialia 8 (2012) 1648–1658

Contents lists available at SciVerse ScienceDirect

Acta Biomaterialia journal homepage: www.elsevier.com/locate/actabiomat

Prediction of permeability of regular scaffolds for skeletal tissue engineering: A combined computational and experimental study S. Truscello a,c, G. Kerckhofs b,c, S. Van Bael a,c, G. Pyka b,c, J. Schrooten b,c, H. Van Oosterwyck a,c,⇑ a

Division of Biomechanics and Engineering Design, Katholieke Universiteit Leuven, Celestijnenlaan 300c PB 2419, 3001 Leuven, Belgium Department of Metallurgy and Materials Engineering, Katholieke Universiteit Leuven, Kasteelpark Arenberg 44 PB 2450, 3001 Leuven, Belgium c Prometheus, Division of Skeletal Tissue Engineering Leuven, Katholieke Universiteit Leuven, O&N I Herestraat 49 PB 813, 3000 Leuven, Belgium b

a r t i c l e

i n f o

Article history: Received 5 July 2011 Received in revised form 6 December 2011 Accepted 12 December 2011 Available online 16 December 2011 Keywords: Bone tissue engineering Computational fluid dynamics Scaffold Laser manufacturing

a b s t r a c t Scaffold permeability is a key parameter combining geometrical features such as pore shape, size and interconnectivity, porosity and specific surface area. It can influence the success of bone tissue engineering scaffolds, by affecting oxygen and nutrient transport, cell seeding efficiency, in vitro threedimensional (3D) cell culture and, ultimately, the amount of bone formation. An accurate and efficient prediction of scaffold permeability would be highly useful as part of a scaffold design process. The aim of this study was (i) to determine the accuracy of computational fluid dynamics (CFD) models for prediction of the permeability coefficient of three different regular Ti6Al4V scaffolds (each having a different porosity) by comparison with experimentally measured values and (ii) to verify the validity of the semi-empirical Kozeny equation to calculate the permeability analytically. To do so, five CFD geometrical models per scaffold porosity were built, based on different geometrical inputs: either based on the scaffold’s computer-aided design (CAD) or derived from 3D microfocus X-ray computed tomography (microCT) data of the additive manufactured (AM) scaffolds. For the latter the influence of the spatial image resolution and the image analysis algorithm used to determine the scaffold’s architectural features on the predicted permeability was analysed. CFD models based on high-resolution micro-CT images could predict the permeability coefficients of the studied scaffolds: depending on scaffold porosity and image analysis algorithm, relative differences between measured and predicted permeability values were between 2% and 27%. Finally, the analytical Kozeny equation was found to be valid. A linear correlation between the ratio U3/Ss2 and the permeability coefficient k was found for the predicted (by means of CFD) as well as measured values (relative difference of 16.4% between respective Kozeny coefficients), thus resulting in accurate and efficient calculation of the permeability of regular AM scaffolds. Ó 2011 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved.

1. Introduction Bone tissue engineering scaffolds are porous materials which assist the healing of large or non-healing bone defects when combined with osteoprogenitor cells and appropriate growth factors [1,2]. A porous scaffold must behave as a carrier for cells and molecules during in vitro culture, and it also has to support the external loads after in vivo implantation [3]. Thus, the scaffold architecture is a crucial factor in fulfilling this combination of mechanical and biological requirements. The mechanical properties (i.e. scaffold stiffness) and the geometrical parameters such as pore shape and size, pore interconnectivity and specific surface ⇑ Corresponding author at: Division of Biomechanics and Engineering Design, Katholieke Universiteit Leuven, Celestijnenlaan 300c PB 2419, 3001 Leuven, Belgium. Tel.: +32 16 327067; fax: +32 16 327994. E-mail address: [email protected] (H. Van Oosterwyck).

area have been shown to influence the success of bone scaffolds [4–6]. The scaffold permeability is a key parameter that best represents all the aforementioned geometrical features [7]. It affects the way in which nutrients (such as glucose) and oxygen disperse through the porous scaffold [8], it influences cell seeding efficiency [9], scaffold degradation [10], three-dimensional (3D) in vitro cell culture [10] and ultimately in vivo bone formation [11]. On the one hand, attaining an optimal value for the permeability is not straightforward, as it depends on the desired biological outcome (e.g. a high-permeability value was found to improve in vivo bone formation [3], whereas a lower value was found to enhance cell seeding efficiency [9]). On the other hand, quantifying and controlling the permeability is crucial to understanding its effect on tissue regeneration. The intrinsic permeability coefficient k of an open porous structure is a measure of the ability of a fluid medium to flow through it and, in the simplest formulation, it can be determined

1742-7061/$ - see front matter Ó 2011 Acta Materialia Inc. Published by Elsevier Ltd. All rights reserved. doi:10.1016/j.actbio.2011.12.021

1649

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

by measuring the amount of fluid (or gas) flowing through the porous material in a given time under an applied external pressure, as captured by Darcy’s law [13].

k ¼ ðQ lLÞ=ðADPÞ

ð1Þ

where Q is the volumetric flow rate, DP is the pressure drop, A and L are the cross-sectional area and the length of the porous material, respectively, and l is the dynamic viscosity of the fluid. It is assumed that Darcy’s Law is valid if the Reynolds number (using the mean pores diameter) is below a limit somewhere between 1 and 10 [14]. Several studies measured the permeability coefficient of porous scaffolds or tissues using perfusion systems, either by applying a constant flow rate with a pump and measuring the pressure drop [15] or by applying a hydrostatic pressure head over the sample and measuring the mass or volumetric flow rate with a weighing scale [16,17]. When studies are limited to those modelling fluid flow and permeability in regular scaffolds, either analytical formulae [18] or homogenization techniques [19–21] are used to determine the permeability tensor. In the first case, the semi-empirical Kozeny equation concisely expresses the relation between the (isotropic) permeability k, the porosity U and the specific surface area Ss as follows:

1 /3 k¼ cK S2s

! ð2Þ

with CK the empirical Kozeny constant. This relation has been frequently extended to estimate the permeability coefficient for different porous media (Carman–Kozeny equation [22] and adapted versions [23]) given its wide range of applicability to all types of soil [24]. The main limitation of this equation is the fact that the determination of the CK coefficient strongly depends on the pore geometry and therefore cannot be estimated a priori. Another limitation of the Kozeny equation is that it is intended for isotropic porous media. Indeed, CK is a scalar value rather than a tensor quantity. Applying the Kozeny equation to capture the flow through an anisotropic medium in a given direction will only yield information on a single component (on the diagonal) of the permeability tensor. As an alternative to these analytical approaches, fluid flow modelling on a simplified geometry of the scaffold (i.e. periodic unit cell) was used to compute the permeability coefficients based on the

(a)

De

(b)

average Stokes flow velocity (neglecting any inertial effects, i.e. for low Reynolds numbers) and calculated in response to an applied pressure gradient. The geometrical features of the unit cell used for building up the model (pores and strut size) can be obtained either from computer-aided design (CAD) of the scaffold [25] or from microfocus X-ray computed tomography (micro-CT) image analysis [12,26]. In the first case, care is required, since the designed scaffold architecture can differ from the manufactured one, depending on the manufacturing technique. This, in turn, can strongly influence the flow field [27]. In the second case, although micro-CT can be considered the gold standard for non-destructive assessment of the scaffold architecture [28–30], dimensional values are dependent on the spatial image resolution and the image analysis algorithms used to extract them. Again, this can bias the modeling results. The objectives of this study are twofold. First, the accuracy of computational fluid dynamics (CFD) models for the prediction of the scaffold permeability is determined. This will be done for three different designed, additive manufactured (AM) open porous Ti6Al4V scaffolds, with the same architecture but with different pore sizes and beam thicknesses. CFD models will be created based on two different geometrical inputs for the internal architecture: either based on the scaffold’s CAD design (further referred to as CAD-based modelling) or derived from micro-CT data of the AM scaffolds (further referred to as micro-CT-based modelling). For the latter, the influence of the spatial image resolution and the image analysis algorithms used to extract the pore and strut dimensions of the scaffold architecture on the predicted permeability will be analysed. The predicted permeability values will be compared with experimental measurements of permeability. Secondly, the validity of the semi-empirical Kozeny equation for analytically calculating the permeability will be verified, and its application within the context of a scaffold design process will be discussed.

2. Materials and methods 2.1. Scaffold design, production and surface treatment The scaffolds used in this study were open porous titanium alloy (Ti6Al4V) structures produced by a non-commercial selective

(c)

Di Top view Beam thickness

Beam length Pore size

H

Side view

Overall dimensions Height H [mm] Internal diameter Di [mm] External diameter De [mm] Fig. 1. (a) CAD of the Ti6Al4V scaffold. (b) Top view and side view of a cross section of the scaffold. (c) Repetitive unit cell.

30 6 12

1650

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

laser melting (SLM) machine [31]. The designed scaffolds had a height H of 30 mm, a diameter Di of 6 mm and a 3 mm thick, fully dense Ti6Al4V annulus (De 12 mm) surrounding the porous structure (Fig. 1). The latter was designed in order to laterally confine the scaffold and to force the fluid to flow through its pores. Three different scaffolds were designed based on the same unit cell (Fig. 1c), but with a different pore size and beam thickness and therefore resulting in three different porosities, further referred to as low (LP), medium (MP) and high (HP) (Table 1).

Table 1 Geometrical features of the scaffold CAD design (see also Fig. 1c); pore size, beam thickness, beam length and porosity values are shown for the LP, MP and HP scaffolds. Scaffold CAD design

Low porosity (LP)

Pore size (mm) Beam thickness (mm) Beam length (mm) Porosity (%)

Medium porosity (MP)

High porosity (HP)

0.75 0.40

0.80 0.30

0.90 0.25

1.15 56.7

1.10 70.5

1.15 79.6

Fig. 2. Representative SEM image of an SLM produced Ti6Al4V scaffold beam surface before (left) and after (right) CHE treatment.

LP

(A)

MP

HP

45 HR

Frequency [%]

40 35 30 25 20 15 10 5 0 0

200

400

600

800

1000

Beam thickness [µm]

(B) 45

270 LR

40

240 210

30

180

25

Grey levels

Frequency [%]

35

20 15 10

150 120 90 60

5

30

0

0 0

200

400

600

Beam thickness [µm]

800

1000

0

100

200

300

400

500

600

Pixels number

Fig. 3. (A) Beam thickness distribution for the LP, MP and HP scaffold resulting from the HR (top) and LR (bottom) micro-CT scan. The average values were used as input for the SF models. (B) Longitudinal micro-CT section of the scaffold of the MP specimen (top) and grey levels (bottom) along the white line drawn on the image. The pore size and the beam thickness were calculated manually by counting the number of pixels (and multiplying by the pixel size) with a grey level either below (for pore size) or above (for beam thickness) a certain threshold.

1651

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

After manufacture, the scaffolds were subjected to a surface treatment, i.e. chemical etching (CHE), in order to reduce the roughness of the strut surface caused by non-fully melted attached powder grains inherent to the manufacturing process (Fig. 2). Briefly, the scaffolds were submerged twice in a stirred solution composed of hydrofluoric acid (HF) and distilled water (0.5 ml HF and 50 ml distilled H2O). The influence of the CHE on the scaffolds’ morphology was tested with several treatment durations (2, 4, 6, 8, 10 and 12 min). Scanning electron microscopy (SEM) based analysis revealed that 10 min of CHE removed all attached grains, after which the scaffolds were cleaned with water and ethanol.

used to estimate the beam thickness distribution. The manufactured beam length was assumed to be equal to the CAD designed beam length (Table 1) and, consequently, the pore size was calculated as the difference between the beam length and the average beam thickness [31]. For the SP method, the pore size and the beam thickness were determined manually based on the grey level of longitudinal micro-CT sections of the scaffold by means of CTan (Fig. 3B). Two scaffold longitudinal micro-CT sections (see side view in Fig. 1b) were analysed for each specimen; five unit cells were selected in different locations of the longitudinal micro-CT section, and one measurement per unit cell was performed.

2.2. Micro-CT-based image analysis of the scaffold geometry

2.3. Permeability analysis

Micro-CT scanning of the manufactured scaffolds, after CHE, was performed using a Philips HOMX 161 X-ray system with AEA Tomahawk CT software [32]. Two different spatial image resolutions were tested: a low-resolution (LR) scan and a high-resolution (HR) scan with an isotropic voxel size of 27.6 and 12.5 lm, respectively. Morphological parameters were determined, after binarization, using commercially available image analysis CTan software (Skyscan NV, Kontich, Belgium). The threshold value was determined based on a grey value plot of a single strut. Instead of having sharp edges, a transition in grey values due to the partial volume effect is present at the edges of the struts. The middle of the s-curve of the grey values along a line perpendicular to the strut edge was chosen as the threshold. Two different methods for the extraction of the pore size and beam thickness were developed, further referred to as the sphere fitting (SF) and side plane (SP) method, respectively. In the SF method, the average beam thickness was obtained from the beam thickness distribution, using the model-independent technique from Hildebrand and Rüegsegger [33] (Fig. 3A). Briefly, spheres were fitted in the solid space of the porous structure, and the diameters of the fitted sphere were

The permeability k for the three scaffold porosities was quantified experimentally and computationally (CFD modelling). In all cases the value of k was calculated using Darcy’s law and derived from the slope of the linear plot of the fluid flow rate as a function of the pressure drop across the scaffold. Finally, it was verified to what extent computational and experimental permeability values follow the Kozeny equation (Eq. (2)), for which scaffold porosity and specific surface area were determined.

(A)

2.3.1. Experimental set-up A set-up was built to measure the scaffold permeability (Fig. 4). Water was used as the fluid and was pumped from an open reservoir to the chamber using a peristaltic pump (IPC-N, Ismatec SA) at two different values of flow rates: 6 and 8 ml min1 (resolution of 0.01 ml min1) corresponding to inlet velocities of 3.5 and 4.7 mm s1, respectively. In order to assess the validity of Darcy’s law, the Reynolds number (Re) was calculated as follows:

Re ¼

qmDpores l

ð3Þ

(B)

Piezometric sensor 10 mm

10 mm

Peristaltic pump

C) DAMPENER

Permeability chamber

Titanium scaffold

Temperature sensor

Reservoir

Fig. 4. (A) Design and (B) image of the permeability chamber. (C) Schematic of the perfusion set-up used for the permeability experiments consisting of a reservoir, temperature sensor, peristaltic pump, dampener, permeability chamber and piezometric sensor connected by tubing. The arrow indicates the flow direction.

1652

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

where q and l are the density and dynamic viscosity of the water, respectively, at 21 °C, v is the mean inlet velocity, and Dpores is the mean pore diameter. For both the flow rates, Re was found to be <10, thus proving the validity of Darcy’s law. A flow dampener (Cole-Parmer Instrument Company Ltd., UK) was connected to the set-up just after the pump to reduce the pulsation intrinsic to the pump and to allow for a continuous flow. An in-house developed piezometric sensor with a resolution of 0.5 mm of water (corresponding to 4.9 Pa) was positioned across the permeability chamber in order to measure the pressure drop across it continuously. The fluid leaving the chamber was re-circulated to the reservoir in a closed loop. A temperature sensor was used to monitor the temperature inside the reservoir, to be able to compensate for any temperature changes, which will in turn affect the viscosity value and, consequently, the permeability coefficient. Two repeated measurements for each scaffold specimen and for each flow rate were carried out. Four measurements (two per flow rate) were conducted for the empty chamber (without any scaffold specimen) in order to determine the pressure drop, resulting from the housing and the connectors at the inlet and the outlet. The scaffold permeability kscaffold was then calculated as follows:

kscaffold ¼

QLl QLl ¼ AðDPscaffold Þ AðDPtotal  DPchamber Þ

ð4Þ

where DPscaffold is the pressure drop across the scaffold specimen, DPtotal is the total pressure drop measured with the scaffold specimen inside the chamber, and DPchamber is that measured for the empty chamber. 2.3.2. Fluid flow modelling Five different 3D geometrical models per scaffold porosity (15 models in total) were built and meshed in the ACIS-based solid modeller Gambit 2.2 (Fluent, Lebanon, NH, USA) as input for the

Inlet velocity=0.1 mm/s

Symmetry condition

No slip condition

CFD analysis. The difference between the models is the way in which the pore size and beam thickness were quantified: they were based either on the CAD design (CAD-based model) or on the micro-CT data of the manufactured scaffold specimens (micro-CT-based). For the latter, a further distinction was made, depending on the spatial image resolution and the method used to extract the dimensional values (see Section 2.2). Hence four different micro-CT-based models were built based on four different geometrical inputs: (i) low-resolution scan with sphere fitting algorithm (LR-SF); (ii) low-resolution scan with side plane method (LR-SP); (iii) high-resolution scan with sphere fitting algorithm (HR-SF); and (iv) high-resolution scan with side plane method (HR-SP). To predict the pressure and velocity fields inside the scaffolds, the commercially available finite volume code Fluent 6.3 (Fluent, Lebanon, NH, USA) was used to set up and solve the fluid dynamic problem. Owing to the symmetry and the repeatability of the scaffold geometry, the domain was limited to one geometrical unit cell in width and two in length (Fig. 5). The edge length of the tetrahedral volumes was equal to 30 lm for all the models. The culture medium was regarded as an incompressible and homogeneous Newtonian fluid with the properties of water (a viscosity of 103 Pa s and a density of 103 kg m3 at a temperature of 21 °C). A constant velocity of 0.1 mm s1 was assigned to the inlet, corresponding to a flow rate of 0.2 ml min1 (Reynolds number Re <1), thus in the range of applicability of Darcy’s law. A zero gauge pressure was applied to the outlet, symmetry conditions to the lateral surfaces and no-slip conditions to the walls of the scaffold. Steady-state Navier–Stokes equations were used to describe the flow problem. The pressure drop across the two unit cells was obtained for all the models and used to calculate the permeability coefficient. 2.3.3. Determination of porosity and specific surface area The porosity values of the manufactured scaffolds were obtained in three different ways: (i) based on micro-CT image analysis and calculated as the percentage of void voxels (determined by thresholding) relative to the total number of voxels; (ii) based on CFD geometrical models and calculated as the ratio between the fluid volume and the total volume; and (iii) experimentally, measured using the Archimedes principle [34]. For the latter, three scaffold specimens per porosity were used, having the same internal architecture as the specimens for the permeability measurements (same pore and beam dimensions) and a similar surface treatment but without the Ti6Al4V annulus. Similarly, the specific surface area was also determined using micro-CT image analysis and the CFD models. In the former case, the specific surface area was calculated as the ratio between the sum of the areas of the pixels at the solid/void interface and the total number of voxels. For the latter, it was calculated as the sum of the area of the elements of the fluid/solid interface (corresponding to the surface wall, grey surface in Fig. 5) divided by the total volume.

3. Results

Zero gauge pressure Fig. 5. Geometry representation of the modelled fluid domain and boundary conditions. The white arrows represent the flow direction. A constant velocity of 0.1 mm s1 was assigned to the inlet (blue surface), zero gauge pressure to the outlet (green surface), symmetry condition on the lateral wall (pink surfaces) and no slip condition to the wall (grey surfaces).

The manufacturing process systematically generated scaffolds with beam thicknesses larger than the designed ones as measured by micro-CT (Table 2). The values of the pore size and beam thickness for the three scaffold porosities (LP, MP and HP) were dependent on the spatial image resolution (LR vs HR) and on the measurement method used (SF vs SP). The biggest difference in pore size between CAD and manufactured scaffolds was found for the low-porosity scaffolds using the LR-SF method, and equalled 22%.

1653

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

Table 2 Dimensional values of pores and beams, as designed (CAD) and as measured by micro-CT image analysis using the four different methods (SP and SF) and spatial image resolutions (LR and HR). For each scaffold porosity, this gives rise to five different CFD models. CFD geometrical models Acronym

Computer aided design method CAD

Low porosity (LP) Medium porosity (MP) High porosity (HP)

Low resolution scansphere fitting algorithm

Low resolution scan-side plane method

LR-SF

LR-SP

Pore size (lm)

Beam thickness (lm)

Pore size (lm)

Beam thickness (lm)

Pore size (lm)

Beam thickness (lm)

Pore size (lm)

Beam thickness (lm)

750

400

583

567

661 ± 36

382 ± 16

654

496

649 ± 18

486 ± 25

800

300

712

388

734 ± 27

343 ± 28

766

384

726 ± 25

365 ± 33

900

250

829

321

814 ± 34

288 ± 28

817

333

795 ± 42

343 ± 38

R

2

=1

Pressure drop Δ P [Pa]

350 300 250

2 R = 0.8795

200

2 R = 0.9546

150 100

2 R = 0.8262

50 0 0

HR-SP

Beam thickness (lm)

Empty chamber LP scaffold MP scaffold HP scaffold

400

HR-SF

High resolution scan-side plane method

Pore size (lm)

The relationship between the experimentally applied flow rate and the measured pressure drop for the three scaffold porosities

450

High resolution scansphere fitting algorithm

2

4

6

8

10

Flow rate Q [ml/min] Fig. 6. Correlation between applied water flow rate (Q) and measured total pressure (DPtotal) difference over the LP, MP and HP scaffolds and the empty chamber (two measurements per flow rate).

(A)

and for the empty chamber is depicted in Fig. 6. As the contribution to the total pressure drop from the components of the set-up (inlet connectors and inlet–outlet channels of the permeability chamber) could not be neglected, the pressure drop over the chamber was taken into account for the measured pressure drop over the three scaffold types. The permeability coefficient kscaffold (mean ± standard deviation) was calculated from Eq. (3) and equalled 0.52  1010 ± 3.9  1012 m2, 1.69  109 ± 1.78  1010 m2 and 3.61  109 ± 8.5  1010 m2, respectively, for the LP, MP and HP scaffold. For the different models (five models per porosity, see also Section 2.3.2), the velocity and pressure distribution were analysed. The pressure at the inlet was calculated as the mean value of the pressure distribution in a cross section perpendicular to the flow direction (y-direction) for x = 0, corresponding to the scaffold inlet (Fig. 7). The pressure value at the outlet was equal to zero for all the models, according to the assigned boundary condition. The resulting value of the pressure drop across the two unit cells and the assigned velocity at the inlet were used to calculate the permeability coefficient, according to Eq. (1). Computed and measured permeabilities can be found in Fig. 8. The permeability values computed from the CAD-based models were substantially higher than those from the micro-CT-based models (relative difference between 43% and 83%) and also than the ones from the experimental tests (relative difference between 131% and 244%). The relative

(B)

Fig. 7. Representative pressure distribution at inlet and outlet cross-sections used to calculate (A) the pressure drop and (B) the velocity distribution in the modelled region of interest for the HR-SP model. The data are shown for the high-porosity scaffold for an inlet velocity of 0.1 mm s1.

1654

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

Permeability coefficient [x 10 -9 m 2 ]

9

CAD based LR-SF model

8

LR-SP model HR-SF model

7

HR-SP model Experimental test

6 5 4 3 2 1 0

LP

MP

HP

Fig. 8. Permeability coefficient as predicted by CFD analysis and as measured experimentally (four measurements per scaffold porosity) for the HP, MP and LP scaffolds. The CFD results correspond to the five models per porosity in Table 2.

CAD based 4

SF-LR model

y = 0,935x

SF-HR model

y = 0,7419x

3,5

SP-LR model

y = 1,0789x

SP-HR model

K CFD [x10-9 m 2 ]

3 2,5 y = 0,3884x

2 y = 0,7613x

1,5 1 0,5 0

0

1

2

3

4

5

6

7

8

9

k experimental [x10-9 m 2 ] Fig. 9. Correlation between the predicted and measured permeability values and linear regression equations.

difference between measured and computed values was strongly reduced when a micro-CT-based model of the manufactured scaffold was applied: between 12% and 43% for the LR-SF models and between 26% and 117% for the LR-SP models, respectively. Relative differences were further reduced using HR scans, resulting in low values between 2% and 27% for HR-SF models and between 13% and 15% for HR-SP models, respectively. Furthermore, the predictive capability of the models was evaluated through a plot relating the computational prediction of scaffold permeability to the experimentally measured values (Fig. 9). Both the HR models were found to give the best prediction of the experimental values (linear regression with a slope equal to 1.08 and 0.93, respectively, for the HR-SP and HR-SF models). The porosity and specific surface values resulting from the CFD models and experimental analyses are shown in Fig. 10. Comparing the experimental porosity values extracted from the micro-CT image analysis with those from the Archimedes principle, a maximum difference of 7% was found both for the LR and HR scans. Comparing porosity values from the Archimedes principle with those from CFD models, the largest differences were found for the LR-SF models (between 1% and 44%). For the other models, differences were between 1% and 28%. Fig. 10B shows good agreement between the specific surface area values obtained from the CFD models and those calculated

experimentally from micro-CT analysis: relative differences between 0% and 26% and between 0% and 21% were found on comparing the CFD model data with the LR and HR scan data, respectively. A linear relation was found between the permeability and the ratio U3/Ss2 for the three scaffold porosities, with U the porosity and Ss the specific surface area (see Eq. (2)). This was the case for the measured as well as the computed (by means of CFD) permeability (R2 = 0.997 for both cases). For the computed values, the HR-SP model was used, as it predicted the best the experimental data (relative difference averaged over the three scaffold porosities equalled 12%). The Kozeny coefficient, being the inverse of the slope of the linear regression curve, was calculated to be 4.8 and 4 for the experimental and computed values, respectively, thus representing a relative difference of 16.4% (Fig. 11). Finally, for all 15 CFD geometrical models, the computed permeability values were plotted as a function of the computed U3/ Ss2, again showing good correlation (R2 = 0.96, see Fig. 12). 4. Discussion The measured permeability values showed very good agreement with the permeability values as computed by means of both

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

1655

Fig. 10. (A) Porosity and (B) specific surface area as determined by the CFD models (geometrical models CAD, LR-SF, LR-SP, HR-SF, HR-SP; see Table 2) and by experimental analysis for HP, MP and LP scaffolds. For the experimental analysis, porosity and specific surface area were calculated based on micro-CT images of the scaffolds (LR and HR scans) and porosity also based on the Archimedes principle.

K

Fig. 11. Correlation between the measured and computed (HR-SP model) permeability as a function of the ratio U3/Ss2; values of porosity U and specific surface area Ss were derived either from micro-CT data (HR scan; for measured permeability) or from a CFD model (HR-SP model; for computed permeability). A table with the values of the Kozeny coefficient (inverse of the slope of the linear correlation curve) for both approaches and their relative difference is shown at the top.

the HR-SP and the HR-SF models. For the former, relative differences of 13%, 15% and 13% were found for the HP, MP and LP scaffolds, respectively, while for the latter relative differences of 8%, 27% and 2% were calculated. The predictive capability of the HRSP and HR-SF model can be further appreciated from the plot of the computed vs the experimental values of permeability (Fig. 9), with a perfect prediction resulting in a linear regression with a slope equal to 1 (and intercept equal to 0). A linear regression with a slope equal to 1.07 and 0.93 was found for the HR-SP and HR-SF models, respectively, thus demonstrating the good predictive capability of these two models. Possible sources of discrepancy between the measured and predicted values can be due to: (i) an insufficient voxel resolution for

the micro-CT data to resolve the surface roughness (see Fig. 2), therefore leading to an underestimation of the specific surface area (its square being inversely proportional to the permeability coefficient); (ii) the viscous effects at the lateral surfaces of the scaffolds confinement which are not taken into account in the models (symmetry boundary conditions applied on the lateral walls), but which can play a role during the experiments, causing an increase in pressure drop and, consequently, a decrease in the measured permeability value; and (iii) the pressure drops occurring at the inlet and outlet regions of the scaffold in the experimental set-up, due to the restriction of the cross-sectional area, which were not considered in the model, and which may lead to an overestimation of the permeability value by means of the CFD predictions. A linear correlation between the ratio U3/Ss2 and the permeability coefficient k was found for the modelling as well as for the experimental approach (Fig. 11), demonstrating the validity of the Kozeny equation (Eq. (2)). The inverse of the slope of the linear correlation curve (Ck constant) is a characteristic of the internal geometry of a porous medium [35]. The CK constant was found to be 4.8 and 4 for the experimental tests and the HR-SP model, respectively (relative difference of 16.4%). Based on the value of this constant, the permeability coefficient of SLM scaffolds, with the same internal porous architecture as the ones studied here, can be calculated analytically (Eq. (2)). Furthermore, the fact that all permeability values, computed by means of all 15 geometrical CFD models, fitted one single linear curve (Fig. 12), demonstrated that differences between these computed values can be explained by the differences in terms of the ratio U3/Ss2 between the corresponding geometrical models. This also means that defining a proper geometrical model should be done with care and must include the porosity and the specific surface area of the actual manufactured scaffolds, hence taking into account production-related constraint with regard to the original scaffold design. For both the porosity values, determined from HR and LR microCT scans and from the Archimedes method, as well as for specific surface area values, obtained from HR and LR micro-CT scans, there

1656

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

Fig. 12. Computed scaffold permeability as a function of the computed ratio U3/Ss2 for the 15 CFD geometrical models (CAD, LR-SF, LR-SP, HR-SF, HR-SP; see Table 2).

were significant differences in function of the method used (Fig. 10). There are several explanations for these differences. First, there is the spatial resolution (voxel size) of micro-CT images, and hence the dimensions of the smallest details that can be captured.

This means that very refined geometrical details may not be captured by the micro-CT images (i.e. details <25 lm cannot be discriminated when a voxel size of 12.5 lm is used), while they are still detected by Archimedes principle [34]. Secondly, the grey

k

k

K

Fig. 13. Flow chart of a scaffold design process for permeability.

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

value threshold for image segmentation was determined manually based on the grey value histogram, but was kept constant, independent of the spatial image resolution, which might also introduce errors. The porosity values based on LR scans were systematically higher than those based on HR scans. This can be attributed mainly to the partial volume effect, causing blurring of the edges between pores and beams. As the grey value threshold was equal for the two scan resolutions, the partial volume effect is expected to be larger for the LR than for the HR scan, as shown in different studies [36,37]. The threshold and the spatial resolution selection of the microCT images may have affected the values of the specific surface area for the same reasons as mentioned above. In particular, the specific surface area values determined from HR scans systematically resulted in higher values than those from LR scans (difference between 10% and 18%). This can be explained by the fact that more geometrical details are captured in the HR images, leading to larger total surface area, while the total volume remains constant. Similar findings were reported in other studies [36,38]. Differences in terms of dimensional values were evident between the different CFD models (Table 2), in turn leading to differences in terms of porosity, specific surface area (Fig. 10) and permeability (Fig. 8). The most substantial dimensional differences were found between the CAD-based models, on the one hand, and the micro-CT-based models on the other. The designed pore size was systematically larger than the pore size of the manufactured specimens (inherent to the SLM technique (see also [31]), leading to a substantial overestimation of the porosity and permeability for the CAD-based models. The differences between the different micro-CT-based models were in general (much) smaller than those between the CAD-based and micro-CT-based models. For the micro-CT-based models, two findings strongly emerged from the comparison between predicted and experimental permeability values. First, the LR scan models led to larger differences with respect to the experimental values (on average 44% difference for the LR scan models, compared with 13% difference for the HR scan models). This was found for both methods (either SP or SF) used to determine the pore size and the beam thickness. This result strengthens the importance of selecting a proper spatial image resolution for the extraction of the scaffold dimensional values, as discussed previously with respect to the porosity and the specific surface area. Secondly, when using HR scans, the differences between the SP and the SF methods was not substantial, and both methods gave rise to models that could accurately predict the permeability values. As the SF method allows for a much more automated determination of the beam thickness and is applicable to any scaffold architecture, it is preferred over the SP method, which makes use of a priori knowledge on the internal porous architecture, but involves manual measurements, making it more labour intensive and time consuming. Therefore, one may prefer to use the SF method for determining the scaffold dimensional values. Finally, the fact that the Kozeny equation is valid for this unit cell geometry leads to a very effective design process for scaffold permeability, as illustrated in Fig. 13. For a given target permeability value, the corresponding target ratio U3/Ss2 can be calculated analytically. For regular scaffolds, this ratio will be a function of one or more architectural parameters (e.g. only pore size, as for this study), in turn leading to target values of these parameters. Knowledge on the systematic offset between manufactured and designed scaffold architecture (such as pore size) finally allows for the establishment of the scaffold CAD design. Although the validity of this approach was demonstrated for a specific unit cell architecture, the entire procedure can be repeated and validated for any unit cell design. This was investigated by performing CFD analyses for three additional unit cell geometries, each for two porosities. CK values

1657

were found ranging between 1.6 and 3.2 (simulation details not included). The advantage of regular scaffolds for tissue engineering applications, previously demonstrated for the prediction of oxygen transport inside scaffolds [39], is further extended here to the assessment of scaffold permeability. The validity of the Kozeny equation and the comparison with experimental values were performed by means of 15 CFD models within a reasonable amount of (computer) time, owing to the fact that only two unit cells along the length and one unit cell along the width of the scaffold were modelled. Thanks to the periodicity of the architecture, this is sufficient to fully represent the flow field for the entire scaffold. For irregular scaffolds, a much larger region of interest should be modelled to be representative of the entire scaffold, leading to much larger computational costs [40]. Furthermore, the approach proposed in Fig. 13 will not apply to irregular scaffolds, owing to the heterogeneity and the lack of control of scaffold architecture. In addition, it is useful to compare the measured and predicted permeability values with those of bone tissue. Kohles et al. [16] measured permeability values for healthy cancellous bone, harvested from mature bovine femur, ranging from 2.33  1010 m2 to 4.65  1010 m2 as a function of the orientation. Permeability values of 8.05  109 m2 and 2.76  109 m2 were predicted for bone specimens harvested from the human vertebral body and the human proximal femur, respectively [41]. The discrepancy between values is due to the structural anisotropy of trabecular bone and the different trabecular microarchitecture in the various anatomical sites. The values of intrinsic permeability of the regular scaffold, which were the subject of the present study (range 5.2  1010 m2 to 3.61  109 m2), are of the same order of magnitude as those of cancellous bone.

5. Conclusion CFD models, based on high-resolution micro-CT images proved to be accurate for the prediction of the permeability of regular AM Ti6AL4V scaffolds, regardless of the method used to derive the geometrical models. For the SP method, relative differences of 13%, 15% and 13% were found for the high-porosity, medium-porosity and low-porosity scaffolds, respectively, while for the SF method these were 8%, 27% and 2%, respectively, in comparison with permeability measurements. It is concluded that high-resolution micro-CT data provided reliable geometrical data to be used as input for CFD analysis for the prediction of scaffold permeability. For the unit cell geometry studied, the validity of the Kozeny equation for analytically calculating the permeability was demonstrated. The method developed can be extended to other unit cell designs, thus leading to an effective scaffold design process.

Disclosures No competing financial interests exist.

Acknowledgements The authors thank Nele Meeuwsen for technical assistance and K.U.Leuven IDO project 05/009 (QuEST) for financial support. G. Kerckhofs was funded by the agency for Innovation by Science and Technology in Flanders (IWT/OZM/080436 and IWT/OZM/ 090655). G. Pyka was funded by the Research Foundation – Flanders (FWO, G.0618.10). This work is part of Prometheus, the Leuven Research and Development Division of Skeletal Tissue Engineering of K.U.Leuven (www.kuleuven.be/prometheus).

1658

S. Truscello et al. / Acta Biomaterialia 8 (2012) 1648–1658

Appendix A. Figures with essential colour discrimination Certain figures in this article, particularly Figs. 1, 3–5, 7–11, and 13 are difficult to interpret in black and white. The full colour images can be found in the on-line version, at doi:10.1016/ j.actbio.2011.12.021. Appendix B. Supplementary data Supplementary data associated with this article can be found, in the online version, at doi:10.1016/j.actbio.2011.12.021. References [1] Cancedda R, Giannoni P, Mastrogiacomo M. A tissue engineering approach to bone repair in large animal models and in clinical practice. Biomaterials 2007;28:4240–50. [2] Porter JR, Ruckh TT, Popat KC. Bone tissue engineering: a review in bone biomimetics and drug delivery strategies. Biotechnol Prog 2009;25:1539. [3] Hutmacher DW, Schantz JT, Lam CXF, Tan KC, Lim TC. State of the art and future directions of scaffold-based bone engineering from a biomaterials perspective. J Tissue Eng Regen Med 2007;1:245–60. [4] Akay G, Birch MA, Bokharib MA. Microcellular polyHIPE polymer supports osteoblast growth and bone formation in vitro. Biomaterials 2004;25:3991–4000. [5] Dennis JE, Haynesworth SE, Young RG, Caplan AI. Osteogenesis in marrowderived mesenchymal cell porous ceramic composites transplanted subcutaneously: effect of fibronectin and laminin on cell retention and rate of osteogenic expression. Cell Transplant 1992;1:23–32. [6] O’Brien FJ, Harley BA, Yannas IV, Gibson LJ. The effect of pore size on cell adhesion in collagen-GAG scaffolds. Biomaterials 2005;26(4):433–41. [7] Li S, De Wijn JR, Li J, Layrolle P, De Groot K. Macroporous biphasic calcium phosphate scaffold with high permeability/porosity ratio. Tissue Eng 2003;9:535–48. [8] Botchwey EA, Dupree MA, Pollack SR, Levine EM, Laurencin CT. Tissue engineered bone: measurement of nutrient transport in three-dimensional matrices. J Biomed Mater Res 2003;67A:357–67. [9] Impens S, Chen Y, Mullens S, Luyten F, Schrooten J. Controlled cell-seeding methodologies: a first step toward clinically relevant bone tissue engineering strategies. Tissue Eng C 2010;16:1575–83. [10] Jeong CG, Hollister SJ. Mechanical, permeability, and degradation properties of 3D designed poly(1, 8 octanediol-co-citrate) scaffolds for soft tissue engineering. J Biomed Mater Res B Appl Biomat 2010;93:141–9. [11] Hui PW, Leung PC, Sher A. Fluid conductance of cancellous bone graft as a predictor for graft–host interface healing. J Biomech 1996;29(1):123–32. [12] Singh R, Lee PD, Lindley TC, Dashwood RJ, Ferrie E, Imwinkelried T. Characterization of the structure and permeability of titanium foams for spinal fusion devices. Acta Biomater 2009;5:477–87. [13] Darcy H. Les Fontaines Publiques de la Ville de Dijon. Paris: Delmont; 1856. [14] Arora KR. Soil mechanics and foundation engineering. (geotechnical 7th engineering) edition. Delhi: Standard Publishers Distributors; 2009. [15] Ochoa I, Sanz-Herrera JA, García-Aznar JM, Doblaré M, Yunos DM, Boccaccini AR. Ceramic TiO2-foams: characterisation of a potential scaffold. J Eur Ceram Soc 2003;24:661–8. [16] Kohles SS, Roberts JB, Upton ML, Wilson CG, Bonassar LJ, Schlichting AL. Direct perfusion measurements of cancellous bone anisotropic permeability. J Biomech 2001;34:1197–202. [17] Jeong CG, Zhang H, Hollister SJ. Three-dimensional poly(1, 8-octanediol–cocitrate) scaffold pore shape and permeability effects on sub-cutaneous in vivo chondrogenesis using primary chondrocytes. Acta Biomater 2001;7:505–14.

[18] Innocentini MDM, Salvini VR, Macedo A, Pandolfelli VC. Prediction of ceramic foams permeability using Ergun’s equation. Mater Res 1999;2:283–9. [19] Sanz-Herrera JA, Garcia-Aznar JM, Doblare M. On scaffold designing for bone regeneration: a computational multiscale approach. Acta Biomater 2009;5:219–29. [20] Hollister SJ, Lin CY. Computational design of tissue engineering scaffolds. Comput Method Appl Mech Eng 2007;196:2991–8. [21] Coelho PG, Fernandes PR, Rodrigues HC. Multiscale modeling of bone tissue with surface and permeability control. J Biomech 2010;44:321–9. [22] Carman PC. Flow of gases through porous media. New York: Academic Press; 1956. [23] Yu B. Analysis of flow in fractal porous media. Appl Mech Rev 2008;61:5081–100. [24] Taylor DW. Fundamentals of soil mechanics. New York: John Wiley; 1948. [25] Lee KW, Wang S, Lu L, Jabbari E, Currier BL, Yaszemski MJ. Fabrication and characterization of poly(propylene fumarate) scaffolds with controlled pore structures using 3-dimensional printing and injection molding. Tissue Eng 2006;12:2801–11. [26] Van Cleynenbreugel T, Schrooten J, Van Oosterwyck H, Vander Sloten J. MicroCT-based screening of biomechanical and structural properties of bone tissue engineering scaffolds. Med Bio Eng Comput 2006;44:517–25. [27] De Boodt S, Truscello S, Ozcan SE, Leroy T, Van Oosterwyck H, Berckmans D, et al. Bimodular flow characterisation in tissue engineering scaffolds using computational fluid dynamics and particle imaging velocimetry. Tissue Eng C 2010;16:1553–64. [28] Lin ASP, Barrows TH, Cartmell SH, Guldberg RE. Microarchitectural and mechanical characterization of oriented porous polymer scaffolds. Biomaterials 2003;24:481–9. [29] Ho ST, Hutmacher DW. A comparison of micro CT with other techniques used in the characterization of scaffolds. Biomaterials 2006;27:1362–76. [30] Jonesa JR, Poologasundarampillaia G, Atwooda RC, Bernardb D, Leea PD. Nondestructive quantitative 3D analysis for the optimisation of tissue scaffolds. Biomaterials 2007;28:1404–13. [31] Van Bael S, Kerckhofs G, Moesen M, Pyka G, Schrooten J, Kruth J-P. MicroCTbased improvement of geometrical and mechanical controllability of selective laser melted Ti6Al4V porous structures. Material Sci Eng 2011;528:7423–31. [32] Kerckhofs G, Schrooten J, Van Cleynenbreugel T, Lomov S, Wevers M. Validation of X-ray microfocus computed tomography as an imaging tool for porous structures. Rev Sci Instrum 2008;79:013711-1–1-9. [33] Hildebrand T, Rüegsegger P. A new method for the model-independent assessment of thickness in three-dimensional images. J Microsc Oxford 1997;185:67–75. [34] Ding M, Odgaard A, Hvid I. Accuracy of cancellous bone volume fraction measured by micro-ct scanning. J Biomech 1999;32:323–6. [35] Sahimi AM. Flow and transport in porous media and fractured rock. Boston: VCH; 1995. [36] Muller R, Koller B, Hildebrand T, Laib A, Gianolini S, Ruegsegger P. Resolution dependency of microstructural properties of cancellous bone based on threedimensional mu-tomography. Technol Health Care 1996;4:113–9. [37] Peyrin F, Salome M, Cloetens P, Laval-Jeantet AM, Ritman E, Ruegsegger P. Micro-CT examinations of trabecular bone samples at different resolutions: 14, 7 and 2 micron level. Technol Health Care 1996;6:391–401. [38] Cooper DML, Matyas JR, Katzenberg MA, Hallgrimsson B. Comparison of Microcomputed Tomographic and Microradiographic Measurements of Cortical Bone Porosity. Calcif Tissue Int 2004;74:437–47. [39] Truscello S, Schrooten J, Van Oosterwyck H. A computational tool for the upscaling of regular scaffolds during in vitro perfusion culture. Tissue Eng Part C Method 2011;17:619–30. [40] Maes F, Van Ransbeeck P, Van Oosterwyck H, Verdonck P. Modeling fluid flow through irregular scaffolds for perfusion bioreactors. Biotechnol Bioeng 2009;103:621–30. [41] Nauman EA, Fong KE, Keaveny TM. Dependence of intertrabecular permeability on flow direction and anatomic site. Annal Biomed Eng 1999;27:517–24.