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