Multiphysics Computational Model Creation
Eight experimental treatments were simulated (static control omitted). The 2 perfusion alone cases (C-P+ and C-P++) were simulated using only Fluent (Workbench 19.0, ANSYS, Inc.), a CFD solver. The remaining six cases were solved using a multiphysics solver, which utilized Fluent, Transient Structural (FE solver), and System Coupling [Fluid-Structure Interface (FSI)] (Workbench 19.0, ANSYS, Inc.). To establish the CFD and multiphysics models, the solid domain was imported to Transient Structural, where the equations of motion was solved in matrix form:
\(\left[M\right]\left\{\ddot{x}\right\}+\left[C\right]\left\{\dot{x}\right\}+\left[K\right]\left\{x\right\}\)=\(\left\{F\left(t\right)\right\}\)
where \(\left[M\right]\left\{\ddot{x}\right\}\)represents the inertia forces,\(\left[C\right]\left\{\dot{x}\right\}\) represents the damping forces, \(\left[K\right]\left\{x\right\}\)represents elastic forces, and\(\left\{F\left(t\right)\right\}\ \)is the dynamic load vector. The material of the solid domain, which is a mixture of PLGA polymer and HA, was assumed to be isotropic elastic. A Young’s modulus and Poisson’s ratio of 2 GPa and 0.3, respectively, were taken from similar PLGA-based scaffolds (Agrawal & Ray, 2001; Holland, Jolly, Yasin, & Tighe, 1987; Zhao et al., 2015). The material density, 2100 kg/m3, was calculated from PLGA (Blasi, D’Souza, Selmin, & DeLuca, 2005) and HA (Larrañaga, Richard J. Lewis, & Lewis, 2007) with 1:1 weight ratio according to the fabrication protocol (Pathi et al., 2010). To simulate the dynamic compression applied to the top of the scaffold (Fig. 1A), a transient 60 μm displacement (5% bulk strain, C+) or 120 μm displacement (10% bulk strain, C++) was applied at the top surface of the solid domain with a 1 Hz sine wave. Bottom surfaces were modeled as fixed supports. The side walls were constrained to allow only in-plane displacements. Finally, fluid-solid shared interfaces were set to allow data transfer (Fig. 3B).
The fluid domain, which filled with Dulbecco’s Modified Eagle Medium, was assumed to be an incompressible fluid with a constant density (1000 kg/m3) and dynamic viscosity (1.45×10-3 Pa\(\bullet\)s) (Olivares, Marsal, Planell, & Lacroix, 2009). In Fluent, Navier-Stokes equations written as a continuity equation:
\begin{equation} \nabla\bullet\left(\rho\overset{}{u}\right)=0\nonumber \\ \end{equation}
and a conservation of momentum equation:
\begin{equation} \frac{\partial\overset{}{u}}{\partial t}+\left(\overset{}{u}\bullet\nabla\right)\overset{}{u}-\nu\nabla^{2}\overset{}{u}=-\nabla\left(\frac{p}{\rho}\right)\nonumber \\ \end{equation}
were solved, where ρ is the flow density, \(\overset{}{u}\) is the flow velocity, p is the static pressure, and ν is the dynamic viscosity. The Pressure-Based Coupled Algorithm of Fluent was used with a convergence criterion of 10-3 for the velocity residual in each direction and for the continuity residual.
In each CFD model, the lateral side boundary was set as a no-slip wall. Boundary conditions of 100 μm/sec velocity and 0 Pa pressure were defined at the bottom inlet and top outlet, respectively, to simulate the direct perfusion in conditions with low perfusion (0.3 mL/min, P+). Similarly, 200 μm/sec bottom inlet velocity was applied to models with high perfusion (0.6 mL/min, P++). To simulate fluid flow induced by compression (C+P-, C++P-), a separate boundary condition in the CFD models was utilized in which both inlet and outlet were instead set at 0 Pa pressure. In order to couple CFD with deforming solid domains, CFD models simulating treatments with compression were set to have deformable mesh with a remesh algorithm. Finally, fluid-solid interfaces were constant (Fig. 3C).
For the multiphysics models, the FE and CFD solvers were linked together by System Coupling in order to perform FSI simulations. During the simulation, the 1 Hz loading cycle was broken into time steps while the equations in either the FE or CFD solver were solved following a staggered iteration strategy. In each time step, the fluid mesh updated its shape following solid deformation, and subsequently, the CFD simulation was solved using the new mesh. Simultaneously, the remesh algorithm was deployed when necessary to avoid highly skewed elements. Then, Transient Structural solved the FE simulation considering the applied boundary and the fluid forces that were transferred from the fluid-solid interface based on the results from the previous CFD calculation. All simulations in this study were performed on remote Linux clusters [The Massachusetts Green High Performance Computing Center (MGHPCC)] with 16 cores [Intel(R) Xeon(R) CPU E5-2650 v3 @ 2.30GHz with 128GB RAM]. After simulations, FE and CFD solvers were linked to PostCFD (Workbench 15.0, ANSYS) for results exporting and visualization, where velocity vectors, WSS color maps, and solid strain color maps were created.
Model Sensitivity Tests
Multiple simulations were carried out for testing model sensitivities. A mesh sensitivity test was performed to determine the appropriate element size. A total of four different models were generated in which the fluid domains were simulated in CFD using the P+ loading condition (low perfusion,100 μm/sec inlet velocity). The models had an element number range from ~44,000 to ~504,000, with the highest element number model was four times smaller in element sizes (both maximum and minimum element lengths) than the lowest. The resulting WSS distributions of each model were nearly identical (Suppl. Fig. 1A), indicating that our starting sensitivity model had sufficiently fine element sizes, similar to other models of tissue engineered scaffold (Zhao, Melke, Ito, van Rietbergen, & Hofmann, 2019). In fact, coarser models did not pass mesh quality tests or had internal errors. We selected the mesh with the same element sizes as our previous CFD study (minimum element length of 0.02 mm, roughly four times of initial micro-CT voxel size) to maintain consistency (Liu et al., 2018).
A second sensitivity test was performed to determine time step length used in multiphysics models. The C+P+ loading condition was applied in three different simulations with 40, 100 and 200 time steps per loading cycle (1 cycle, 1 second). Distributions of fluid-solid interface WSS of each model at 0.25 s were virtually identical (Suppl. Fig. 1B). Thus, we chose 100 time steps per loading cycle to achieve a balance between data points and simulating time. Next, we assumed the scaffold was sufficiently stiff to prevent the fluid from causing reciprocal solid deformations. To verify this, we conducted two-way FSI coupled models to compare the solid strain distributions to the compression-alone FE simulation. Identical solid strain distributions verifyied that fluid forces did not deform solid structures (Suppl. Fig. 1C). Thus, we utilized one-way coupling simulations, allowing only solid deformation transfers. Finally, to determine the number of loading cycles to simulate, the C+P+ loading condition was used to simulate two full loading cycles (2 sec). WSSs were calculated from the fluid-solid interface on each loading cycle, and median values of these WSS showed identical curves between the first and the second cycle (Suppl. Fig. 1D). Thus, 1 second of loading was simulated.
Results