Geometry-Based Computational Fluid Dynamic Model for Predicting the Biological Behavior of Bone Tissue Engineering Scaffolds

The use of biocompatible and biodegradable porous scaffolds produced via additive manufacturing is one of the most common approaches in tissue engineering. The geometric design of tissue engineering scaffolds (e.g., pore size, pore shape, and pore distribution) has a significant impact on their biological behavior. Fluid flow dynamics are important for understanding blood flow through a porous structure, as they determine the transport of nutrients and oxygen to cells and the flushing of toxic waste. The aim of this study is to investigate the impact of the scaffold architecture, pore size and distribution on its biological performance using Computational Fluid Dynamics (CFD). Different blood flow velocities (BFV) induce wall shear stresses (WSS) on cells. WSS values above 30 mPa are detrimental to their growth. In this study, two scaffold designs were considered: rectangular scaffolds with uniform square pores (300, 350, and 450 µm), and anatomically designed circular scaffolds with a bone-like structure and pore size gradient (476–979 µm). The anatomically designed scaffolds provided the best fluid flow conditions, suggesting a 24.21% improvement in the biological performance compared to the rectangular scaffolds. The numerical observations are aligned with those of previously reported biological studies.


Introduction
Modern tissue engineering strategies rely on the use of additive manufacturing with biocompatible and biodegradable materials, living cells and growth factors to create constructs to restore, repair, and improve the function of damaged tissues and/or organs [1][2][3]. There are two commonly used strategies, cell-laden and scaffold-based approaches [4][5][6]. The cell-laden approach is based on the use of additive manufacturing to selectively print bio-inks (hydrogels containing cells) [7,8]. A range of hydrogels has been investigated, and the fabrication of multi-material and cellular constructs has been reported [9,10]. However, cell-laden constructs have limited mechanical properties, making them more suitable for soft-tissue applications [11,12]. Scaffolds are the most common approach used in bone tissue engineering [4,5].
Scaffolds are 3D biocompatible and biodegradable porous structures that provide the necessary substrate for cell attachment, proliferation, differentiation, extracellular matrix (ECM) formation, and mineralisation [6,7]. They must be highly porous structures

•
Use of different blood flow velocities (BFV) through scaffolds that are used for largesized defects located in bones [37][38][39]. • Newtonian fluid: fluid is assumed to have a constant viscosity to simplify the system. It is noted that WSS value distribution can change when using a non-Newtonian fluid [40,41]. • Surface roughness: The surface was assumed to be smooth based on previous experimental data, as it is dependent on the fabrication method. It also must be noted that including roughness can change how the fluid interacts with the scaffolds [42].

Scaffold Design
Two scaffold configurations were considered for simulations purposes based on previously published experimental work for two scaffold design configurations [43,44]. The first configuration (Case 1) corresponds to rectangular scaffolds designed with a 0/900 lay-down pattern ( Figure 1) and three different pore sizes (Case 1A:300 µm pore size; Case 1B:350 µm pore size; Case 1C:450 µm pore size). The scaffolds were designed with 12 filament layers, 350 µm in diameter, and 330 µm in layer thickness. The second configuration (Case 2) corresponded to anatomically designed scaffolds with six rings ( Figure 2) and a graded pore size (first ring, 476 µm; second ring, 629 µm; third ring, 670 µm; fourth ring, 730 µm; fifth ring, 803 µm; sixth ring, 979 µm). The scaffold porosity was determined using the following equation: where V scaffold is the volume of the scaffold filaments, and V solid is the volume of the solid part. The specific surface area (S), which defines the available area for cell attachment, proliferation, and differentiation and influences permeability (increases with increasing pore size and decreasing S) was calculated as follows: where SA is the surface area and V is the volume. The geometric characteristics of the scaffolds are listed in Table 1.
where Vscaffold is the volume of the scaffold filaments, and Vsolid is the volume of the solid part. The specific surface area (S), which defines the available area for cell attachment, proliferation, and differentiation and influences permeability (increases with increasing pore size and decreasing S) was calculated as follows: where SA is the surface area and V is the volume. The geometric characteristics of the scaffolds are listed in Table 1.

Simulation
CFD simulations were conducted to investigate the fluid flow characteristics of different scaffolds, including permeability, WSS, and flow velocity. In contrast to other CFD studies that have considered water as the fluid, this study used blood as the fluid medium ( Table 2). The fluid was considered incompressible and Newtonian, and the scaffolds were modelled as rigid bodies within the fluid domain [36][37][38]. Moreover, the 3D Navier-

Simulation
CFD simulations were conducted to investigate the fluid flow characteristics of different scaffolds, including permeability, WSS, and flow velocity. In contrast to other CFD studies that have considered water as the fluid, this study used blood as the fluid medium ( Table 2). The fluid was considered incompressible and Newtonian, and the scaffolds were modelled as rigid bodies within the fluid domain [36][37][38]. Moreover, the 3D Navier-Stokes equations were considered to solve the conservation of mass and momentum across the scaffold using the finite volume method. The steady state Navier-Stokes equations are given by [45]: where ρ represents the fluid density (kg/m 3 ), t is the time (s), u is the fluid velocity (m/s), µ is the dynamic viscosity (Pa.s), ∇ represents the del operator, p is the pressure (Pa) and F denotes other forces (i.e., gravity or centrifugal force) acting in the fluid domain. In this study, F is assumed to be zero [46,47]. The intrinsic permeability of the scaffolds, k 0 , (m 2 ), for a range of inlet velocities is given by [48]: where L is the scaffold length (m) and ∆P is the pressure drop (Pa). However, it should be noted that for very high inlet velocities, Darcy's law is inadequate for estimating the permeability of the scaffold [24,[39][40][41]. Typically, when the calculated interstitial Reynolds number (Re) is higher than 8.6, Darcy's classic equation (Equation (5)) is invalid [39,42]. Therefore, to evaluate the validity of Darcy's law, the Reynolds number (Re) was determined as follows [28,41]: where d is the pore size (m). Assuming a laminar flow regime, the WSS can be calculated throughout the domain for each scaffold by considering the normal velocity gradient over a filament [49,50]: where τω represents the wall shear stress (Pa) and n indicates the x-, y-, and z-directions (m). CFD simulations were conducted using the finite volume method (FVM) using the ANSYS CFX solver (ANSYS INC, Washington, PA, USA). The scaffolds were designed using SolidWorks (Dassault Systems, Vélizy-Villacoublay, France) and the CAD models were imported into the ANSYS CFX design modeller, where the filaments were fused to create a single body. A cylindrical enclosure was created around the scaffold to simulate a flow in a bone line environment. The fluid domain was constructed via a Boolean subtraction operation using the scaffold as a tool and a cylinder with a radius of 0.001 m and 0.01 m on either side of the centre of the scaffold.
The inlet was defined at one end of the cylindrical domain, with an outlet on the opposite side. The external region of the flow domain is defined as a rigid-bodied wall (no-slip condition). The domain was meshed using CFX's built-in mesher using the patchconforming method with tetrahedral elements. A mesh independence study was conducted for the three domain refinements, resulting in differences of less than 5%. The total number of elements for the considered scaffolds was 3,441,553 for Case 1A, 4,180,282 for Case 1B, 4,154,355 for Case 1C, and 6,384,572 for Case 2.
To mimic physiological conditions, both isothermal heat transfer at 37 • C and a shear stress transport (SST) turbulence model were considered. The SST model is a two-equation k-w model that improves the prediction of separation in the near-wall regions of the model [45]. Therefore, the model can account for turbulence which is an physiological characteristic of blood flow [51].
In the pre-processor the inlet was defined with medium (5%) turbulence, and a range of inlet blood flow velocities (1-9 mm/s) were used following experimental conditions simulating tibia bone fractures [52]. The outlet was set as a 0-resistance boundary condition. The convergence criteria for the residuals were set to a root-mean-square value of 1e-4 using second order methods. Each simulation was completed on an Intel-Xeon quad-core CPU using the MPI. The results were post-processed using the CFX CFD-post program to extract the velocity, pressure, and WSS values. A histogram of the WSS versus the geometrical percentage was also exported. Figure 3 shows the pressure drop variation as a function of inlet velocity for the different scaffolds. Results were obtained by considering two virtual probes on the top and bottom of the scaffolds and finding the corresponding value from a plane contour passing through the scaffold. As observed, the highest-pressure drop was obtained for the Case 1A and Case 1C scaffolds, followed by Cases 1B and 2. These results can be explained by the combined effect of pore size and architecture of the scaffold (number of filaments).

Results and Discussion
The high pressure drop observed for Case 1A can be explained by the larger interaction between the fluid and the scaffold due to the large number of filaments, high specific area, and the smaller diameter, while for Case 1C the results are due to the low interactions (among the rectangular scaffolds this is the configuration presenting the lowest number of filaments and the lowest specific surface area) and large pore size. The gradient architecture in Case 2 explains the lowest pressure drop. The high pressure drop observed for Case 1A can be explained by the larger interaction between the fluid and the scaffold due to the large number of filaments, high specific area, and the smaller diameter, while for Case 1C the results are due to the low interactions (among the rectangular scaffolds this is the configuration presenting the lowest number of filaments and the lowest specific surface area) and large pore size. The gradient architecture in Case 2 explains the lowest pressure drop.
Using the pressure drop values, it was possible to estimate permeability using Darcy's law, and the results showed that Case 2 scaffolds presented the highest permeability, followed by Cases 1 B, 1A, and 1C ( Figure 4). Moreover, the different in permeability decreases with the increase in blood flow velocity, indicating a larger effect of pore size and geometries at lower velocity compared to higher velocity. Nevertheless, the trend is consistent with different velocities applied. As observed, the anatomically designed scaffolds (Case 2) guaranteed sufficient blood flow through the scaffolds, whereas among the rectangular scaffolds, the best results were achieved with Case 1 B. Using the pressure drop values, it was possible to estimate permeability using Darcy's law, and the results showed that Case 2 scaffolds presented the highest permeability, followed by Cases 1 B, 1A, and 1C ( Figure 4). Moreover, the different in permeability decreases with the increase in blood flow velocity, indicating a larger effect of pore size and geometries at lower velocity compared to higher velocity. Nevertheless, the trend is consistent with different velocities applied. As observed, the anatomically designed scaffolds (Case 2) guaranteed sufficient blood flow through the scaffolds, whereas among the rectangular scaffolds, the best results were achieved with Case 1B. Figures 5-16 show the changes in the fluid metrics (velocity, pressure, and WSS) as a function of the inlet velocity. As observed, for rectangular scaffolds, the pressure increases with increasing pore size (Case 1C presents the highest-pressure values followed by Cases 1A and Case 1B) (Figures 5-7). However, in the case of scaffolds with a gradient of pore sizes (Case 2), regions with higher average pore sizes presented lower pressure values (Figure 8). For all considered scaffolds, results show that an increase in the inlet velocity results in larger values of pressure due to the larger deflection forces occurring when the fluid encounters the scaffolds of higher pore sizes and lower specific surface areas. For the velocity, the results suggest no impact of pore size (Figures 9-12).
For the WSS (Figures 13-16), the maximum value was observed in Case 1B, followed by Cases 1C, 1A, and 2. In Cases 1A and 1B, high WSS values were observed at the top of the scaffold, with a significant drop from the top surface to the inner zone. However, in Case 1C, high WSS values were observed on both the top and side faces of the scaffold. This can be attributed to the scaffold geometry and specific area, suggesting a correlation between a higher WSS and lower specific surface area. As WSS has a major impact on the biological performance of scaffolds, the results suggest that the specific surface area and pore size are key design parameters. Moreover, the maximum pressure was incident on the outer edges of all scaffolds; however, it was relatively lower for Case 2, which explains the overall lower WSS values observed on the anatomically designed scaffolds ( Figure 16). Therefore, to reduce the intensity of pressure on the sides of the scaffolds, Case 2 scaffolds can be used because they present circular geometries and a higher specific area, resulting in a reduction in the pressure at the top surface and overall pressure. Thus, the results suggest that the Case 2 scaffolds present an improved overall biological behavior compared to the more commonly used rectangular scaffolds (Case 1).  Figure 4. The effect of inlet velocity on the permeability. * Statistical evidence (p < 0.05) analysed by one-way ANOVA, and Tukey post hoc test. The * statistical evidence (p < 0.05), **, *** and **** is the one-way analysis of variance (one-way ANOVA) and Tukey's post hoc test with the use of GraphPad Prism software and is used to show the difference between the results. The * is a small difference, while more * are added as the differences between the results increases. Figures 5-16 show the changes in the fluid metrics (velocity, pressure, and WSS) as a function of the inlet velocity. As observed, for rectangular scaffolds, the pressure increases with increasing pore size (Case 1C presents the highest-pressure values followed by Cases 1A and Case 1B) (Figures 5 to 7). However, in the case of scaffolds with a gradient of pore sizes (Case 2), regions with higher average pore sizes presented lower pressure values ( Figure 8). For all considered scaffolds, results show that an increase in the inlet velocity results in larger values of pressure due to the larger deflection forces occurring when the fluid encounters the scaffolds of higher pore sizes and lower specific surface areas. For the velocity, the results suggest no impact of pore size (Figures 9-12).               For the WSS (Figures 13-16), the maximum value was observed in Case 1B, followed by Cases 1C, 1A, and 2. In Cases 1A and 1B, high WSS values were observed at the top of the scaffold, with a significant drop from the top surface to the inner zone. However, in Case 1C, high WSS values were observed on both the top and side faces of the scaffold. This can be attributed to the scaffold geometry and specific area, suggesting a correlation between a higher WSS and lower specific surface area. As WSS has a major impact on the biological performance of scaffolds, the results suggest that the specific surface area and pore size are key design parameters. Moreover, the maximum pressure was incident on the outer edges of all scaffolds; however, it was relatively lower for Case 2, which explains the overall lower WSS values observed on the anatomically designed scaffolds ( Figure 16). Therefore, to reduce the intensity of pressure on the sides of the scaffolds, Case 2 scaffolds can be used because they present circular geometries and a higher specific area, resulting in a reduction in the pressure at the top surface and overall pressure. Thus, the results suggest that the Case 2 scaffolds present an improved overall biological behavior compared to the more commonly used rectangular scaffolds (Case 1).        To further analyse the WSS distribution, the percentage of the total surface area of the scaffolds with WSS values lower than 30 mPa against the inlet velocity was also considered ( Figure 17). This WSS value was used as a reference, as it is the cut-off value for cell proliferation and differentiation. As observed for all scaffolds, for velocities ranging between 1 and 5 mm/s the occurrence of WSS > 30 mPa was very low, between 5 and 7 mm/s the occurrence of WSS > 30 mPa slightly increased, and for velocities ranging between 7 and 9 mm/s the occurrence of WSS > 30 mPa significantly increased. The highest WSS values were observed for Case 1C, followed by Cases 1A, Case 1B and Case 2. The average WSS values are listed in Table 3. In the case of rectangular scaffolds, the results suggest that the highest average cell viability is expected to occur in Case 1B scaffolds (the scaffold percentage of volume with WSS < 30 mPa is 61.64%), followed by Case 1C (61.56%), and Case 1A (59.44%). Among all scaffold architectures, the highest cell viability was expected to occur in Case 2 (85.86%), with the highest WSS values at the sides and the lowest WSS at the centre.  To further analyse the WSS distribution, the percentage of the total surface area of the scaffolds with WSS values lower than 30 mPa against the inlet velocity was also considered ( Figure 17). This WSS value was used as a reference, as it is the cut-off value for cell proliferation and differentiation. As observed for all scaffolds, for velocities ranging between 1 and 5 mm/s the occurrence of WSS > 30 mPa was very low, between 5 and 7 mm/s the occurrence of WSS > 30 mPa slightly increased, and for velocities ranging between 7 and 9 mm/s the occurrence of WSS > 30 mPa significantly increased. The highest WSS values were observed for Case 1C, followed by Cases 1A, Case 1B and Case 2. The average WSS values are listed in Table 3. In the case of rectangular scaffolds, the results suggest that the highest average cell viability is expected to occur in Case 1B scaffolds (the scaffold percentage of volume with WSS < 30 mPa is 61.64%), followed by Case 1C (61.56%), and Case 1A (59.44%). Among all scaffold architectures, the highest cell viability was expected to occur in Case 2 (85.86%), with the highest WSS values at the sides and the lowest WSS at the centre.
These results agree with experimental studies that investigated the biological behavior of similar rectangular scaffolds, showing that the Case 1B scaffolds presented the best results [43,44]. Other experimental studies also agree with the improved biological behavior that anatomically designed scaffolds (Case 2) present compared to rectangular scaffolds (Case 1) [43,44]. The strategy of this study differs from other published works in that it uses a simple design by controlling the pore size and geometry by changing the fibre size and distance. Other studies have used a unit cell approach that may have an irregular geometrical flow, as they are based on metallic AM. However, our findings agree with the literature using similar assumptions that permeability and WSS are geometry-dependent, and that there is no clear trend with pore sizes and geometries.
This strategy is based on previously published experimental studies using extrusionbased additive manufacturing with polymers and polymer composites. The use of anatomically designed scaffolds (Case 2) is a new and promising strategy to replace conventional rectangular scaffolds. Previously published experimental and numerical results obtained in this study can be used to optimize different scaffold geometries, reducing the time and cost of repeated biological experiments. Future perspectives may include the use of graded porous scaffolds and non-Newtonian and transient simulations to further develop numerical strategies for bone tissue engineering scaffolds.  These results agree with experimental studies that investigated the biological behavior of similar rectangular scaffolds, showing that the Case 1B scaffolds presented the best results [43,44]. Other experimental studies also agree with the improved biological behavior that anatomically designed scaffolds (Case 2) present compared to rectangular scaffolds (Case 1) [43,44].

Conclusions
The strategy of this study differs from other published works in that it uses a simple design by controlling the pore size and geometry by changing the fibre size and distance. Other studies have used a unit cell approach that may have an irregular geometrical flow, as they are based on metallic AM. However, our findings agree with the literature using similar assumptions that permeability and WSS are geometry-dependent, and that there is no clear trend with pore sizes and geometries.
This strategy is based on previously published experimental studies using extrusionbased additive manufacturing with polymers and polymer composites. The use of anatomically designed scaffolds (Case 2) is a new and promising strategy to replace conventional rectangular scaffolds. Previously published experimental and numerical results obtained in this study can be used to optimize different scaffold geometries, reducing the time and cost of repeated biological experiments. Future perspectives may include the use of graded porous scaffolds and non-Newtonian and transient simulations to further develop numerical strategies for bone tissue engineering scaffolds.

Conclusions
Different scaffold geometries have been investigated using Computational Fluid Dynamics (CFD) to understand the effect of pore size and geometry on biological behavior which is usually difficult to obtain experimentally. Among the rectangular scaffolds, Case 1B exhibited the optimum balance between pore size and specific surface area, showing the highest permeability, WSS < 30 mPa at higher velocities, and a lower pressure drop. Anatomically designed scaffolds (Case 2) presented similar velocities but improved WSS distribution, suggesting superior biological behavior. The numerical results were aligned with previously obtained experimental results.
In the future, a simulation model will be developed using 3D scanning data of printed scaffolds to better describe the architecture of the scaffolds, including surface roughness which may impact the results. Additionally, transient simulations, advanced fluid-structure interaction models, and advanced fluid rheology models, including the use of pseudoplastic models for the blood model, were considered.