Wall Shear Stress Analysis and Optimization in Tissue Engineering TPMS Scaffolds

When designing scaffolds for bone tissue engineering (BTE), the wall shear stress (WSS), due to the fluid flow inside the scaffold, is an important factor to consider as it influences the cellular process involved in new tissue formation. The present work analyzed the average WSS in Schwartz diamond (SD) and gyroid (SG) scaffolds with different surface topologies and mesh elements using computational fluid dynamics (CFD) analysis. It was found that scaffold meshes with a smooth surface topology with tetrahedral elements had WSS levels 35% higher than the equivalent scaffold with a non-smooth surface topology with hexahedral elements. The present work also investigated the possibility of implementing the optimization algorithm simulated annealing to aid in the design of BTE scaffolds with a specific average WSS, with the outputs showing that the algorithm was able to reach WSS levels in the vicinity of 5 mPa (physiological range) within the established limit of 100 iterations. This proved the efficacy of combining CFD and optimization methods in the design of BTE scaffolds.


Introduction
In bone tissue engineering (BTE), scaffolds are porous support matrixes designed as an environment to promote cell proliferation, differentiation, and growth [1]. These cellular processes are considerably influenced by several parameters, one of which is the wall shear stress (WSS) that affects the cells inside the scaffold [2][3][4][5]. WSS is caused by the relative movement between the scaffold walls and the fluidic phase inside the scaffold. Studies have shown that different levels of WSS lead to different mechanical signals for the mesenchymal stromal cells, which, in turn, cause differences to the cellular differentiation process [6,7]. This parameter is affected by various factors, with one being the scaffold geometry. In fact, small changes in the scaffold's design parameters (such as its porosity, surface topology, or curvature) considerably alter the average WSS experienced by the cells inside the scaffolds [8,9]. Furthermore, in order to promote the desired conditions to promote bone proliferation, the average WSS experienced by the cells needs to be between 0.1 and 10 mPa [10,11].
When designing scaffolds, the multiple possible inputs for scaffold development result in a large range of possible designs. Therefore, computational simulations are normally implemented to reduce the time and cost associated with conceptualizing a new scaffold. For BTE, computational fluid dynamics (CFD) analysis is essential in understanding how changes to a scaffold's design influence its fluidic properties (such as the WSS) and the underlying cellular processes [12][13][14].
Additionally, computational simulations can also be used to design scaffold geometry with the desired characteristics, through a structural optimization approach. Although there are several different optimization methods, most of them follow the same framework [15]. After defining the design space, the chosen objective function, and the constraints for that specific optimization process, the process starts with an initial step where the initial geometry (initial solution) is set as well as the material properties of the structure. Afterward, the state variables are computed through a proper numerical tool in order to analyze the relevant properties of the initial geometry. If the resulting properties do not reach the objective, then an optimization algorithm is utilized to obtain a new geometry, according to the predefined constraints. This process repeats itself iteratively until a new structure satisfies the objective function. However, in BTE, optimization methods have rarely focused on the fluidic properties of the scaffold [14][15][16], instead only consider the scaffold's mechanical properties (such as their Young's modulus [17]; compressive strength [18] and octahedral shear strain [19][20][21]).
Taking this into consideration, the objective of this work was twofold. Firstly, the difference in average WSS was analyzed for different scaffold topologies, with either smoothed (tetrahedral elements) or non-smoothed (hexahedral elements) wall surface topologies. This was in order to determine whether previously designed numerical models [22] can be used in the study of WSS, or if an additional step to correct their surface topology is required. Secondly, a simulated annealing (SA) optimization strategy was used to design scaffolds for a specific average WSS, in order to investigate the feasibility of using structural optimization methodologies combined with CFD analysis in the design of scaffolds for required fluidic properties.

Scaffold Design
To generate the scaffold's mesh for the CFD analysis a computer code was used, which, given the specific design parameters, created a hexahedral mesh of the desired geometry for the analysis [22]. This mesh consists of a single cubic unit of the fluidic phase of the chosen scaffold design, with forty elements per side, seeing as previous studies have shown that 40 elements per side are the minimum in order to obtain an appropriate scaffold geometry [23]. For this study, the chosen scaffold geometries were the triply periodic minimum surfaces (TPMS) gyroid (SG) and Schwartz diamond (SD) designs. This choice was because these scaffolds were found to be the most appropriate TPMS structures for cellular growth given their high permeability and fluid tortuosity that promotes cellscaffold interaction [22]. Additionally, an empty portion was added before and after the scaffold, which represented the empty permeability chamber that allows the fluid flow to stabilize before and after passing the scaffold (Figure 1), mimicking the experimental setup for permeability analysis [24].
Having designed a scaffold mesh with hexahedral elements (Figure 2a), to create the corresponding mesh with tetrahedral elements and a smooth surface, an additional design step was required. To this end, Meshlab ® [25] was used to apply a three-step Laplacian smoothing process followed by a three-step isotropic explicit remeshing process to obtain the scaffold with smooth surfaces (Figure 2b).
The scaffold design parameters considered for this paper were: the porosity of the scaffold and the length of the side of a single basic cubic unit (higher values result in scaffolds with larger pore sizes). In addition to these two characteristics, two other parameters that were analyzed were the length of empty chambers before and after the scaffold and, in the case of the tetrahedral meshes, the world unit parameter from Meshlab's isotropic explicit remeshing process (which directly controls the number of elements in the final mesh). Given that these factors were constant across all scaffolds that were studied, a convergence test was carried out for each one to determine their best value. For these tests, both scaffold geometries were analyzed in regard to their average WSS and permeability with a fluid inlet velocity of 0.0001 m/s, a length of 5 mm for the cubic unit, and a 70% porosity. geometries were analyzed in regard to their average WSS and permeability with a fluid inle velocity of 0.0001 m/s, a length of 5 mm for the cubic unit, and a 70% porosity.

Surface Topology CFD Analysis
Having defined the optimal values for the cubic unit length and the world unit pa rameters, the next step was to define the scaffold designs to be compared. For each of th two TPMS geometries, three different levels of porosity were considered (60%, 70%, and 80%) for a total of six different scaffold designs. A hexahedral and a tetrahedral mesh wer then created for each design and their average WSS was studied using CFD analysis wit the previously mentioned fluid inlet velocity of 0.0001 m/s and cubic unit length of 5 mm For this study, the computational simulations were conducted using the commercia  geometries were analyzed in regard to their average WSS and permeability with a fluid inlet velocity of 0.0001 m/s, a length of 5 mm for the cubic unit, and a 70% porosity.

Surface Topology CFD Analysis
Having defined the optimal values for the cubic unit length and the world unit parameters, the next step was to define the scaffold designs to be compared. For each of the two TPMS geometries, three different levels of porosity were considered (60%, 70%, and 80%) for a total of six different scaffold designs. A hexahedral and a tetrahedral mesh were then created for each design and their average WSS was studied using CFD analysis with the previously mentioned fluid inlet velocity of 0.0001 m/s and cubic unit length of 5 mm. For this study, the computational simulations were conducted using the commercial

Surface Topology CFD Analysis
Having defined the optimal values for the cubic unit length and the world unit parameters, the next step was to define the scaffold designs to be compared. For each of the two TPMS geometries, three different levels of porosity were considered (60%, 70%, and 80%) for a total of six different scaffold designs. A hexahedral and a tetrahedral mesh were then created for each design and their average WSS was studied using CFD analysis with the previously mentioned fluid inlet velocity of 0.0001 m/s and cubic unit length of 5 mm. For this study, the computational simulations were conducted using the commercial software Fluent ® Ansys ® (Ansys Inc., Canonsburg, PA, USA), which has already proven to be a useful tool for analyzing the fluidic properties of scaffolds [26,27]. In terms of the parameters for the Fluent solver, the fluid chosen for the simulations was water, in accordance with previous studies [22,28], with a density of 1000 kg/m 3 and dynamic viscosity of 0.001 Pa.s. The fluid was also assumed to pass through the scaffold in a steady state laminar flow.
Three different boundary conditions need to be defined, these being the velocity inlet, the pressure outlet, and the wall boundary. For these simulations, it was considered that the fluid was traveling in the y direction through the scaffold. This meant that an exterior wall on the xz plane was the velocity inlet surface, while the opposite wall was the pressure outlet with 0 Pa, so that the pressure drop of the scaffold would be equal to the pressure at the inlet (since the pressure drop is the pressure at the inlet minus the pressure at the outlet). Additionally, interfaces with periodic boundary conditions were implemented in parallel on the remaining exterior walls in order to create an infinite scaffold in the x and z directions.
Finally, after running the simulations, the software CFD-post from Ansys was used to calculate the average WSS on the scaffold walls and, if needed, the pressure difference between the inlet and outlet of each scaffold. Afterward, the scaffold permeability could be calculated using the pressure drop and Darcy's law (Equation (1)) [29,30]: where K is the permeability expressed in m 2 ; ∆P is the pressure drop expressed in Pa; L is the length of the scaffold expressed in m; v is the inlet velocity of the fluid and µ is the dynamic viscosity of the water which is 0.001 Pa*s.

Structural Optimization Process
For this work, the created optimization algorithm was focused on the scaffold's average WSS on both the SD and SG geometries with tetrahedral elements. More specifically, the algorithm will attempt to design a scaffold geometry with an average WSS that promotes cellular proliferation and growth.  [11] discussed how the WSS of structures meant for bone growth needed to be between 0.1 and 10 mPa. Given that the present method looks at the average WSS, a target was set at 5 mPa for the optimization program, as an average target value between the two edges of the accepted physiological levels. Additionally, seeing as the average WSS of the scaffold was dependent on the fluid velocity passing through it (unlike permeability, which is only dependent on the scaffold design), a constant inlet velocity of 0.001 m/s (similar to the parameter described in a comparable scaffold perfusion analysis [31]), was selected to allow a comparison between all scaffold geometries.
To solve the established scaffold optimization challenge, the metaheuristic optimization computational algorithm simulated annealing (SA) was chosen. This is a metaheuristic optimization approach used to solve optimization problems with a large search space. Additionally, this algorithm allows a worse solution to be accepted in the earlier iterations, which minimized the probability of the program reaching a local minimum and being unable to progress further.
SA requires the definition of an appropriate cooling schedule variable (a); the objective function; the maximum number of iterations the algorithm runs before stopping (K max ); and the upper (U j ) and lower (L j ) boundaries of each parameter. In this work, a cooling schedule variable of 0.9 was chosen, alongside a maximum of 100 iterations. The objective function was defined as the normalized difference between the target (5 mPa) and the real average calculated WSS, and the design variables are the side length of the unit cell of the TPMS structure and the scaffold porosity. These variables are limited to the lower and upper bounds of 1 and 10 mm for the unit length and 60% and 80% for the porosity, respectively. These values were chosen to maintain pore sizes that promote cellular processes.
After defining the variables, the algorithm then calculated an initial value for each parameter (equal to the halfway point between the upper and lower boundary of each parameter), and used the process discussed in the previous section to create and study the average WSS of the scaffold.
The program then used a Matlab ® (Mathworks, Natick, MA, USA) script to calculate the objective function value of the designed scaffold and recorded it as the current solution (Z c ). Subsequently, the algorithm defined an initial temperature (T j ) equal to 1/5 of the current solution. The temperature was cooled after each iteration by multiplying it by the previously established cooling schedule variable (a).
Having the initialization of the SA algorithm, the iterative process started with testing the immediate neighbors of the current solution (using the current parameters and normal distribution). If a parameter went outside the predefined boundaries, then it was replaced by the closest values within the boundaries (for example, if the new porosity parameter was equal to 85%, then the code replaced it with 8', the corresponding maximum allowed value). The program then calculated the objective function value of the new parameters, using the same method as discussed above, and recorded it as the next trial solution (Z n ). If Z n was lower than or equal to Z c , then the code always accepted the new value. Otherwise, the chance of accepting the new value was equal to e ((Z c −Z n )/T j ) . The code then recorded the value of the parameters and objective function value and proceeded with the temperature cooling. Finally, when the algorithm ended, it returned the history of the optimization process. Figure 3 shows the results of the convergence analysis on the pressure drops and average WSS for simulations conducted on both 70% porosity test scaffolds. After obtaining these results, the relative difference between the analysis with the smallest world unit and every other simulated mesh was calculated (Tables 1 and 2). As expected, the results showed that an increase in the number of elements in the simulation led to more accurate results. However, an increase in the number of elements also resulted in longer simulations. Therefore, a compromise had to be made between the accuracy of the results and the computational cost of the simulation. Accordingly, the world unit value 0.15 was decided as the best option seeing as it always presented a considerably low relative error (almost always around or below 1%) with considerably quicker simulations than models with a higher element count. For the fluid flow to stabilize, the chambers before and after the scaffold must have sufficient "empty" length to allow the flow to stabilize before and after the scaffold; otherwise, this could affect the numerical simulation, namely the pressure drop measurement. However, the longer the empty chamber, the higher the computational cost of the simulations, which would in turn considerably affect the computational cost of the optimization algorithm. Therefore, the desirable empty chamber length is the shortest one that  For the fluid flow to stabilize, the chambers before and after the scaffold must have sufficient "empty" length to allow the flow to stabilize before and after the scaffold; otherwise, this could affect the numerical simulation, namely the pressure drop measurement. However, the longer the empty chamber, the higher the computational cost of the simulations, which would in turn considerably affect the computational cost of the optimization algorithm. Therefore, the desirable empty chamber length is the shortest one that still allows the stabilization of the fluid flow. Furthermore, an additional length should be added to the minimum length to account for different scaffold configurations that require a longer chamber to stabilize, such as lower porosity scaffolds.

Scaffold Design Parameters
In order to determine a minimal length, two different scaffolds were designed: one where the length of the empty chambers, before and after the scaffold, was equal to the length of the basic cubic unit and another where the length empty chambers was equal to half the length of the basic cubic unit. As the fluid flow pressure drop measurement and the average scaffold WSS are both independent from the length of the chambers (as long as the flow is stable at both ends of the numerical model), if both CFD simulations returned similar pressure drop and average WSS values, then this meant that the model with the smaller empty chamber was enough to allow the fluid flow to stabilize.
As shown in Table 3, for both the SG and SD scaffolds, the difference between the models with varying values of h was minimal, with their difference always below 0.3%. This indicates how the shorter chamber was sufficient in allowing the flow to stabilize at the edges of the numerical model. On top of that, analyzing the fluid flow streamlines presented in Figure 4 highlights how the fluid flow is completely stable at the edges of the numerical model. However, choosing a shorter chamber could potentially risk the CFD simulations for the scaffolds with the smallest pore sizes whilst not significantly contributing to a lower computational cost. Therefore, an empty chamber length of half the length of the scaffold was chosen for all the CFD analysis.

Surface Topology CFD Analysis
Tables 4 and 5 illustrate how the difference in the of scaffold surface topology (original jagged surfaces with hexahedral elements or smooth surfaces with tetrahedral) influences the average WSS of the six different scaffold designs. The results show that all of the smoothed scaffolds had a higher average WSS than the original scaffolds (average of 35% increase). Additionally, the volumes of the original and smoothed scaffolds were compared to determine if the smoothing process would influence the porosity of the scaffolds (Table 6).

Surface Topology CFD Analysis
Tables 4 and 5 illustrate how the difference in the of scaffold surface topology (original jagged surfaces with hexahedral elements or smooth surfaces with tetrahedral) influences the average WSS of the six different scaffold designs. The results show that all of the smoothed scaffolds had a higher average WSS than the original scaffolds (average of 35% increase). Additionally, the volumes of the original and smoothed scaffolds were compared to determine if the smoothing process would influence the porosity of the scaffolds (Table 6).

Optimization Method
For each scaffold geometry, the optimization algorithm ran multiple times. Examples of one optimization run for the SD and SG geometries are shown in Tables 7 and 8, respectively. The iterations where no new point was considered were not presented in the tables. The convergence of the results is shown in Figure 5.

Optimization Method
For each scaffold geometry, the optimization algorithm ran multiple times. Exam of one optimization run for the SD and SG geometries are shown in Tables 7 and spectively. The iterations where no new point was considered were not presented i tables. The convergence of the results is shown in Figure 5.

Discussion
The comparison between the different surface topologies showed that the WSS levels in the smoothed surfaces increased by an average of 35% when compared to the original Materials 2022, 15, 7375 9 of 11 non-smoothed surfaces. A previous study has shown that in scaffolds with equal geometry but different surface roughness levels, a lower average WSS was calculated for the scaffolds with a higher roughness [32]. This corresponds to what was observed in Table 5, where the original jagged mesh, given its higher surface roughness, resulted in a lower WSS than the smoothed tetrahedral mesh. Furthermore, when fabricating scaffolds using 3D printing techniques, these methods normally create surface topologies closer to the smoothed surfaces than the original jagged ones. Therefore, the smoothed surface was chosen as the basis for the optimization algorithm of the scaffold's average WSS. Finally, the difference between the volume before and after the smoothing is illustrated in Table 6. The difference between the two is minimal, indicating that the smoothing process can be used without it affecting the scaffold's porosity.
In terms of the optimization algorithm, it was able to consistently reach the desired average WSS of 5 mPa taking a maximum of 12 iterations, well within the established 100-iterations limit. However, given the inherent random nature of the SA algorithm, the number of iterations and the final parameters varied between each optimization run. Nevertheless, the outputs remained similar to those reported in Tables 7 and 8. The optimal SD scaffold was always close to the maximum limit of 10 mm for the cubic unit length and a porosity close to 70%; whilst the SG scaffolds had a cubic unit length between 7 and 8.5 mm and porosity between 63% and 73% (higher porosity corresponded with lower cubic unit length).
The resulting scaffold designs from this algorithm were able to refine the scaffold geometry in order to reach the desired average WSS of 5 mPa. This highlights how the optimization algorithm presented in this work can be combined with CFD analysis; additionally, this algorithm can enhance the design of TPMS scaffolds with a specific fluidic property such as a specific average WSS or permeability. The results also demonstrate how the SA method worked as intended, with the code accepting worse solution to the current one at earlier iterations (for example the SG optimization from the second to the fifth iteration as seen in Table 8) and only accepting strictly better solutions at higher iterations (when the temperature variable T j of the algorithm was at its lowest).
Although the SA algorithm proved as useful to optimize and improve one specific parameter, it would not be able to simultaneously improve multiple fluidic parameters, such as both WSS and permeability. In other words, WSS and permeability are inversely correlated, with studies have shown that that increasing the average WSS leads to a decrease in the scaffold's permeability and vice versa [33,34]. To address the limitations of the current algorithm, a multi-objective optimization approach might be implemented in future developments. Instead of focusing on a single objective function, multi-objective optimization could analyze each parameter individually and would give a set of possible solutions, organized as a Pareto front, but with a higher computational cost. This method would result in multiple possible solutions, which would need to be analyzed in terms of their suitability for BTE.

Data Availability Statement:
The generated models can be made available upon reasonable request to the corresponding author.