InsightsintoInterstitial Flow, Shear Stress,and Mass TransportEffects on ECM heterogeneity in Bioreactor-Cultivated Engineered Cartilage Hydrogels
Tony Chen, Mark Buckley, Itai Cohen, Lawrence Bonassar, Hani A. Awad
Supplemental Information
1. Biphasic theory solution of unidirectional flow over a permeable hydrated mixture
Reproduced from Hou et al. (1989)
Hou et al. (1989) reported a theoretical solution to the classical problem of Poiseuille flow over a porous mixture (Figure A1) using the biphasic theory. Here we reproduce the solution in dimensionless form for reference.
Nomenclature:
Volume fraction of the solid phase in an immiscible biphasic mixture
Volume fraction of the fluid phase in an immiscible biphasic mixture
Flow velocity
Free fluid viscosity (in the flow channel)
Apparent fluid viscosity in the biphasic layer
Darcian permeability
kHydraulic permeabilityk
Interstitial drag coefficient in the biphasic layer, which is inversely proportional to the Darcian permeability
Pressure drop driving the flow in the channel
Height of flow channel
Thickness of biphasic layer
Figure S1.Schematic representation of the problem of Poiseuille flow over a biphasic layer.
Reproduced with modification from Hou et al. (1989).
In setting up the problem, the flow of an incompressible Newtonian fluid (e.g. viscous)is assumed to be driven by a unidirectional pressure gradient along the flow channel and over a biphasic layer, which in turn is assumed to be an immiscible mixture of a linearly elastic, and isotropic solid phase and an incompressible, viscous fluid phase such that . The boundary conditions at the interface between the free fluid and the biphasic layer surface were derived by considering the balance laws for mass, linear momentum, and energy.
The force balance for the fluid in the channel and biphasic layer yields the following governing equations, respectively:
subject to the “pseudo no-slip” and the momentum “jump” boundary conditions at the biphasic layer in contact with the free fluid
and the no-slip boundary condition at the rigid impermeable walls
The governing equations (S1 and S2) can be written in nondimensional form as
Where
is the nondimensional fluid flux
is a nondimensional (thickness) coordinate
is the ratio of the flow channel height to the biphasic layer thickness
is the ratio of viscous drag of the outside fluid to the interstitial drag
of the fluid within the porous biphasic layer
isa weighted viscosity ratio
The pseudo no-slip condition and momentum jump boundary condition at the channel-biphasic layer interface (), respectively, reduce to the following nondimensional form:
The no slip boundary conditions at the rigid impermeable walls
The solution to the nondimensional governing equations subject to the aforementioned boundary conditions can be expressedin a compact form as
Therefore, the interstitial flux depends upon , , and , whose effects have been discussed in detail by Hou et al. (1989). In our setup, it is reasonable to assume that is constant since the free fluid and interstitial fluid are identical and thus have the same viscosity value. The experimentally measured flow velocity values were fitted using least square nonlinear regression (fminsearch function in MATLAB) to the equation , where is the piecewise biphasic solution expressed in equation S11, to determine the values of and . While the nominal dimensions of the chamber and hydrogel were known, we chose to include as an optimization variable to eliminate any experimental uncertainty in determining the z-coordinate of the hydrogel surface. The drag coefficient and the apparent permeability where then determined from the optimized value of as discussed in the manuscript.
2. Computational fluid dynamics modeling
Since the theoretical solution described in section 1 was derived for Poiseuille flow over an infinite slab of a porous material which does not account for the finite geometry of the bioreactor system and the end conditions that might give rise to two-directional, reversible flow, a 2-D computational fluid dynamics (CFD) model of the parallel plate bioreactor system was created using the geometry and dimensionsof the scaled-down bioreactor used in the FRAP-based flow visualization experiments (Table 1). CFD modeling was performed using COMSOL Multiphysics v4.0a (Comsol, Stockholm, Sweden). The bioreactor was discretized into 3regionscomprising the flow media channel, TEC agarose hydrogel, and the anchoring cell-free agarose hydrogel (Figure 2a). To optimize the computational time without sacrificing the accuracy of the numerical solution, the mesh within the TEC hydrogel region, whichcomprised111,864 triangular elements, was refined to be biased towards the interfaces with the free fluid channel and the anchoring agarose layer (Figure S2). The remainder of the geometry was meshed using 105,064 uniformly distributed triangular elements. “No-slip” boundary conditions were enforcedat the walls. The viscous fluid in the media flow channel was modeled using the Navier-Stokes equations, while the hydrogel regions were modeled using the Brinkman’s equation.The TEC hydrogel permeability was set to the value determined for each flow rate (Table 2). The permeability of the cell-free agarose hydrogel was scaled as a function of the agarose concentration as described by Gu et al (2003). The CFD model was solved using the PARDISO static solver. Convergence of the model was confirmed using mesh refinement.
The fluid velocity profilesand flow streamlines (Figure S3) were examined at the midline of the bioreactor and locations 1 mm, either upstream or downstream, from the midline. These measurements were done to verify that the flow is mostly unidirectional at the bioreactor midline, where the FRAP measurements were taken. The model was also used to compute the shear stress profile at midline through the depth of the bioreactor (Figure 5).
Figure S2.(a) The bioreactorcomputational grid showing the different regions comprising the CFD model including the flow channel, the TEC hydrogel, and the anchoring agarose hydrogel layer. (b) Magnification of a region at the interface between the flow channel and the TEC hydrogel showing the biased grid density with increased elements at the interface.
e
Figure S3.Flow velocity profiles within the flow chamber and the TEC hydrogels were estimated from a CFD model of the scaled-down parallel plate bioreactorfor the case of 0.14 ml/min at the midline of the bioreactor and locations 1 mm, either upstream or downstream, from the midline. Panels a andb depict the predicted parallel (y-velocity) and perpendicular (z-velocity) components of the flow velocity, respectively. Note the difference in the scale of the two velocity components, which suggests that the assumption that the flow is mostly unidirectional is valid. Panels c and d show the depth dependent interstitial flow y-velocity and z-velocityprofiles, respectively, within the TEC hydrogel region. Panel e depicts representative streamline plots showing the predicted flow paths within the different bioreactor regions.