|Home | About | Journals | Submit | Contact Us | Français|
Intravascular gas embolism may occur with decompression in space flight, as well as during cardiac and vascular surgery. Intravascular bubbles may be deposited into any end organ, such as the heart or the brain. Surface interactions between the bubble and the endothelial cells lining the vasculature result in serious impairment of blood flow and can lead to heart attack, stroke, or even death. To develop effective therapeutic strategies, there is a need for understanding the dynamics of bubble motion through blood and its interaction with the vessel wall through which it moves. Toward this goal, we numerically investigate the axisymmetric motion of a bubble moving through a vertical circular tube in a shear-thinning generalized power-law fluid, using a front-tracking method. The formulation is characterized by the inlet Reynolds number, capillary number, Weber number, and Froude number. The flow dynamics and the associated wall shear stresses are documented for a combination of two different inlet flow conditions (inlet Reynolds numbers) and three different effective bubble radii (ratio of the undeformed bubble radii to the tube radii). The results of the non-Newtonian model are then compared with that of the model assuming a Newtonian blood viscosity. Specifically, for an almost occluding bubble (effective bubble radius = 0.9), the wall shear stress and the bubble residence time are compared for both Newtonian and non-Newtonian cases. Results show that at low shear rates, for a given pressure gradient the residence time for a non-Newtonian flow is higher than that for a Newtonian flow.
The formation of gas emboli, or bubbles, in blood vessels occurs during deep sea diving, during decompression in long-duration space-flight missions and extravehicular activity, in cardiac surgery, in endoscopy, and in tissue biopsy. Such intravascular gas emboli may deposit into organs, such as the heart or the brain, and as a result cause permanent injury. The distribution of gaseous emboli in the blood vessels is influenced not only by their differential density from blood (buoyancy) but also by the prevailing blood flow dynamics.1 The simplest possible representation of an intravascular gas embolus is a pressure-driven motion of a liquid through a tube containing a freely suspended bubble. In this context, several researchers have investigated (both experimentally and numerically) bubble motion in an arterial blood vessel. 2–5 When the bubble nearly occludes the tube, which is often the case in gas embolism, both its mobility and shape are affected by the presence of confining wall. The bulk phase in gas embolism problems is blood, which is a rheologically complex fluid that exhibits shear-thinning, non-Newtonian characteristics at low shear rate conditions.6 Finite fluid inertia of the blood furthermore complicates the dynamics. These factors necessitate the use of a full numerical solution for describing the bulk and dispersed phase dynamics. The phases are typically characterized by large jumps in properties (density, viscosity) in addition to the presence of high and possibly variable surface tension forces.
The shear-thinning rheology of blood can be adequately described as a function of shear rate by some of the commonly used mathematical/empirical models, namely, Casson, power-law, Carreau and its derivatives Cross and Carreau–Yasuda, Walburn–Schneck, and the generalized power-law models.7–9 Although the Walburn–Schneck and power-law models predict the behavior of blood well at low shear rates, they predict decreasing viscosities at high shear rates and hence fail to capture the Newtonian nature of blood at such high shear rates.8 However, the Walburn–Schneck model has the blood hematocrit (volume fraction of red blood cells) as a model parameter and can thus explicitly account for its possible variations. For example, the latter aspect can take into account the hematocrit variations due to the Fahraeus–Lindquist effects.10 The Casson and the Carreau–Yasuda models capture the experimental behavior of blood well at both high and low shear rates. Furthermore, a modified Casson model can also be used to explicitly account for the hematocrit in the blood.11 The generalized power-law model encompasses the power-law model at low shear rates, a New-tonian model at high shear rates, and for a given hematocrit value, the Casson model as a special case.9 Hence, for this study, blood rheology is adequately described by the generalized power-law model for various shear rates at a hematocrit of 0.45, which is the normal value for humans. The Fahraeus–Lindquist effect tends to reduce the value of tube hematocrit (average hematocrit across a given cross section). However, for tubes in the size range 20–500 µm, there is a formation of a plasma or cell-depleted layer adjacent to the cylindrical vessel wall where the hematocrit is nearly zero.6,10,12 In such vessels, this results in the formation of two distinct flow regions: a clear Newtonian plasma layer adjacent to the cylindrical wall and a central non-Newtonian core consisting of plasma and red blood cells whose hematocrit (core hematocrit) value is slightly higher than the tube hematocrit. The value of the core hematocrit is assumed to be 0.45 in our simulations. For vessels of size larger than 500 µm, the thickness of the plasma layer is assumed to be negligible and the tube hematocrit is taken to be the same as core hematocrit.
For multiphase motion, the results of many studies on transport phenomena associated with the motion of a drop or bubble involving non-Newtonian fluids have been discussed in the monograph of Chhabra.13 Some recent theoretical and numerical works include the creeping motion of shear-thinning power-law fluids past a fluid sphere,14,15 the numerical simulations of bubble motion freely rising in various shear-thinning fluids obeying a power-law model,16 and a Carreau–Yasuda model.17 However, pressure-driven flows appear to have not been considered in any of the non-Newtonian flow studies. With respect to drop motion in pressure-driven flows, many theoretical and numerical studies exist appropriate for low Reynolds number regime of a Newtonian fluid in a cylindrical tube with the drop size comparable to the tube diameter (e.g., Refs. 18–20). For this report, we have numerically simulated the axisymmetric motion of an initially spherical gas bubble moving through a pressure-driven flow of blood modeled as a shear-thinning generalized power-law fluid in a vertical circular channel. To the best of our knowledge, there seem to be no numerical studies that have focused on the behavior of nearly occluding bubbles in a shear-thinning medium where the effects of both pressure-driven flow and gravity are taken into account.
In the context of modeling of multiphase flows, several numerical schemes are available. These methods include front-tracking/immersed boundary,21 level set,22 phase-field, 23,24 volume-of-fluid,25–27 coupled level-set and volume of fluid,28 immersed interface,29 ghost-fluid methods,30,31 and moving mesh interface tracking.32 In all the methods except those listed in Refs. 30–32, flow discontinuities are smoothed and the surface tension force is distributed over a thin layer near the interface to become a volume force. The Navier–Stokes equation is then solved on a fixed Eulerian mesh.
Our numerical study is based on a front-tracking method coupled with a level contour reconstruction procedure33 for periodic redistribution of the front points and reconstruction of a new front. In the front-tracking method, an explicit background mesh of interconnected marker points is used to represent the interface. Such a Lagrangian representation of the interface allows an accurate calculation of the surface tension forces without the direct computation of the interface curvature.
We consider the axisymmetric motion of a bubble in a vertical tube of circular cross section of radius R and height L as shown in Figure 1. L is chosen in such a way that the bubble attains its terminal velocity and shape after release and continues to remain so for at least one bubble diameter. On the basis of numerical experimentation, it was found that L = 8R was sufficient to capture the steady-state dynamics of the bubble for the parametric range considered. The initial radius of the bubble is r0 and the bubble phase fluid is assumed incompressible with density ρg and viscosity μg. The bulk fluid is assumed incompressible with density ρl and viscosity μl. The equations of motion for an isothermal, incompressible two-phase flow problem under study can be expressed by a single fluid continuum model as follows:
In the preceding equations, u is the fluid velocity, p is the dynamic pressure (total pressure minus the hydrostatic head), g is gravitational acceleration, s is the arc length measure on the interface, κ is the curvature of the interface, σ is the constant surface tension coefficient, S(t) denotes the interface (time dependent), nf denotes the unit normal vector on the interface (pointing into the bulk fluid), xf denotes the position vector on the interface, and δ(x − xf) stands for the delta function that is nonzero only when x = xf . The density and viscosity of the medium can be expressed as follows:
where, H(r,z,t) is the Heaviside (step) function, which is 1 in the bulk phase and 0 inside the bubble phase.
The bulk fluid is modeled using a two-fluid model to account for the presence of a cell-depleted layer or plasma layer near the tube wall, whose thickness, ε, depends on tube diameter.6,10,12 As stated in Section 1, the plasma layer is prominent in tubes of sizes 20–500 µm. For tubes of sizes larger than 500 µm, ε is assumed to be negligible and the fluid is assumed to be non-Newtonian everywhere. The plasma layer viscosity μp = 1.5 cP (plasma viscosity) with hematocrit value is equal to zero. Outside the plasma layer, the fluid is assumed to be shear thinning and is modeled using a generalized power-law model9 with a core hematocrit value of 0.45 (average value of human hematocrit). This value is assumed to be the same for the blood vessel sizes considered in this study. The bulk viscosity for the two-fluid model is given as follows:
Here, λ and n are the consistency coefficient and power-law index for the shear-thinning fluid and are assumed to be functions of the local strain rate · μ∞ is the infinite shear rate viscosity (Newtonian viscosity). The Δμ, Δn, a, b, c, and d are the generalized power-law model parameters. For blood, the following values are assumed9: μ∞ = 0.035 g/cm−s, n∞ = 1.0, Δμ = 0.25, Δn = 0.45, a = 50, b = 3, c = 50, and d = 4. is related to the strain rate tensor (D) and is given as follows:
tr(D2) denotes the trace of the tensor D2. The jump in viscosity across the plasma layer in Eq. (4) is smoothed out over two grid cells in our numerical simulation by using a smoothed Heaviside function. Such smoothing is also more physical because under physiological conditions, the boundary between the core and the plasma layer is not sharp and is defined on only an average sense.
For the bubble phase, properties of air are assumed for ρg and μg.
Inlet (z = 0): A fully developed velocity profile corresponding to the generalized power-law fluid is specified by numerically solving an axisymmetric channel flow problem given by the following equation:
In Eq. (7), uz is the velocity component along the axial direction (z axis) and μl is given by Eqs. (4) and (5). Here, the pressure gradient dp/dz is varied until the desired inlet centerline velocity Umax is obtained.
Outlet (z = H): Outflow boundary condition with zero normal derivatives for the velocities is prescribed:
Symmetry boundary (r = 0): A reflective symmetry boundary condition is prescribed at r = 0,
Here, ur is the radial component of velocity.
Wall boundary (r = R): A no-slip boundary condition is prescribed at the wall r = R, given by
As mentioned earlier, the viscosity jump in Eq. (4) is smoothed across the boundary between the core and plasma layer, and hence explicit boundary conditions (continuity of velocity and shear stress) are not considered across this layer for numerical simulations.
At t = 0, a fully developed velocity profile obtained by numerically solving Eq. (7) is prescribed everywhere in the domain. The bubble is initially positioned a diameter away from the inlet boundary.
Nondimensionalization of the governing equations is based on the following scales: 2R for length, Umax, the maximum inlet velocity for the non-Newtonian fluid (inlet centerline velocity), for velocity, 2R/Umax for time, μ∞ for viscosity, ρl for density, and max for pressure. With this scheme, the following nondimensional parameters that govern the motion are obtained:
Here Re characterizes the ratio of inertial to viscous forces, Ca characterizes the ratio of viscous to surface tension forces, We characterizes the ratio of inertial to surface tension forces, and Fr characterizes the ratio of inertial to gravitational forces. From Eq. (5) it can be noted that as → ∞, μl → μ∞ for the shear-thinning fluid. The steady state velocity of the bubble is represented by U∞ in our subsequent discussions.
The method used in this study is based on a front-tracking procedure. The interface/front is explicitly represented by a set of marker points on which the surface tension forces are evaluated. The Navier–Stokes equation is solved on a fixed Eulerian mesh. In the bubble motion problem, the ratios of the properties of the two fluids across the interface/front can be large, for example, of the order of 1000 for density and 100 for viscosity. Such sharp jumps in properties across the interface of ideally zero thickness complicate the numerical simulation because of instabilities. Also, the integration of the surface tension delta source term in Eq. (2) introduces numerical difficulties if the sharp discontinuities are not resolved. To alleviate these problems, all the discontinuities in fluids’ properties are first smoothed out across a finite thickness about the exact interface (diffused interface), and the thickness of this region is proportional to the mesh size. Smoothing of properties is achieved using a phase indicator function (the Heaviside function) whose value varies smoothly from 1 in the continuous (bulk) phase to 0 in the bubble phase. The smooth variation of the phase indicator function across the interface is obtained by solving a Poisson equation (see Refs. 21 and 34 for more details).
Surface tension forces evaluated on the interface are distributed to the Eulerian cells through a numerical approximation of the delta function and by using a “density-weighted” distribution to avoid numerical instabilities caused by high surface tension forces.
When the interface deforms to a large extent or if the volume of the bubble is less than 0.5% of the original volume, an interface reconstruction is undertaken using a level-contour reconstruction procedure.34 Doing so ensures a smooth distribution of front points as well.
The unsteady Navier–Stokes equations (Eqs. (1) and (2)) are then discretized using a finite difference–based variable-density projection method on a MAC-type (Marker And Cell type) grid described in Refs. 35–37. The velocity, density, and viscosity are all located at cell centers. The lagged pressure p(n−1/2) is located at cell corners, with the superscript n denoting the time level. The time stepping procedure is based on a second-order Crank–Nicholson method. The details of the steps involved in discretization of the governing equations for an axisymmetric case and the corresponding time step restrictions for numerical stability are given in Ref. 28. Also, in all our simulations a grid size of 64 × 512 was used, which resulted in grid-independent solutions.
Difficulties with convergence occur for low shear rate conditions where the invariant tends to zero near the core. Consequently, from Eq. (4), the viscosity tends to high values, resulting in convergence problems. In such cases, a lower bound is fixed for (= min) such that a stable convergent solution is obtained. min was varied over a range for which steady-state solutions exist and wherein the differences in such solutions are negligibly small. From our numerical experimentation, a value of min = 0.2 resulted in stable converged solutions with negligible loss in accuracy.
Results are presented for two widely different inlet Reynolds numbers, Re = 0.06, 40. The corresponding characteristic shear rates, Umax/R, are 10 and 180, respectively. A characteristic shear rate of 10 would correspond to blood behavior in the non-Newtonian regime, whereas a value of 180 would essentially correspond to a Newtonian behavior.7 The values of Re chosen are such that they correspond to flow conditions in blood vessels of different sizes, the smaller one corresponding to that of an arteriole and the larger one to that of a small artery. As stated earlier, the axial pressure gradient, dp/dz (in Eq. (7)), is adjusted to obtain the desired inlet centerline velocity, Umax (and hence the desired Re), for the non-Newtonian fluid. In the absence of a bubble this value of dp/dz would result in a fully developed flow profile as given in Figure 2. For the same dp/dz, the results are then compared with that of a Newtonian fluid of viscosity μ∞ (blood assumed to be Newtonian). For a given Re and tube radius R, the initial radius (r0) of the bubble is varied over a specific range, 0.5 ≤ (r0/R) ≤ 0.9, to study the effects of occlusion.
In surface tension–dominated regimes, which occur when values of Re 1, Ca 1 and for Re > 1, We 1, there results the appearance of spurious vortices near the interface called “parasitic currents” that may destroy the interface with time. These nonphysical vortices are numerical artifacts arising from the discretization of the surface tension forces in diffuse interface methods, such as the front-tracking scheme. The surface tension σ for an air–blood interface is approximately 50 dynes/cm. For Re = 0.06 and σ = 50 dynes/cm, Ca = 6.9 × 10−5, and We = 4.0 × 10−6. For Re = 40 and σ = 50 dynes/cm, Ca = 7.8 × 10−3, and We = 0.32. Under such circumstances, the value of surface tension (σ) may have to be artificially reduced up to a value where the effects of parasitic currents are minimized. From our numerical experiments, for Re = 0.06, a value of σ = 0.05 dynes/cm, and for Re = 40, σ = 25 dynes/cm, were found to be suitable and resulted in negligible parasitic vortices for the entire range of r0/R studied. The resulting values of nondimensional parameters for Re = 0.06 are Ca = 6.9 × 10−2, We = 4.0 × 10−3, Fr = 2.25 × 10−2 and for Re = 40 are Ca = 1.6 × 10−2, We = 0.64, Fr = 1.02. The solution for Re = 0.06, Ca = 0.07 and Re = 40, We = 0.6 may still be considered appropriate to describe the dynamics of an air bubble because of surface tension’s still being the dominant force with low Ca and We for both the Res.
This case corresponds to a low shear rate condition (Umax/R = 10), where blood behaves as a non-Newtonian fluid. The size of the tube (diameter) is taken to be 200 µm. As stated earlier, for this tube size, the size of the plasma layer (ε) is assumed to be nonzero and is taken to be 4 µm.12 The nondimensional value of dp/dz = 1100, and the same value is used for comparison with Newtonian fluid of viscosity, µ∞.
In Figure 2, the inlet velocity profiles for the two-fluid generalized power-law model and the Newtonian model are displayed. For consistency, both the velocities are scaled using the inlet maximum velocity (Umax) of the non-Newtonian fluid. The non-Newtonian model reveals a blunted velocity profile (pluglike profile) at the central core of the tube as opposed to a parabolic profile of the Newtonian fluid. For the same dp/dz, the centerline velocity of the Newtonian model is approximately twice that of the non-Newtonian model. A nearly flat velocity profile for the power-law fluid gives rise to a very low value of strain rate that results in a high value of viscosity (see Eq. (4)). This results in increased viscous dissipation and a reduced velocity for the non-Newtonian model compared with the Newtonian case for a given pressure gradient.
In Figure 3, the steady-state bubble shapes are plotted for Re = 0.6 and aspect ratios of 0.5, 0.7, and 0.9, respectively. The vertical sidewalls of the bounding boxes denote the positions of the actual walls of the vertical tube. The axial coordinate is rescaled such that the position of the front stagnation is taken as the zero reference axial point. For r0/R = 0.5, the bubble remains nearly spherical because of high surface tension forces, whereas for r0/R = 0.7, the bubble profile shows slight elongation at the leading meniscus. For r0/R = 0.9, bubble elongation is pronounced. Such pronounced deformation seen here may be the result of an artificial reduction of surface tension compared to the actual value. The curvature at the leading meniscus is larger than the trailing meniscus because of the applied pressure gradient. The direction of the applied pressure gradient is such that a higher pressure acts at the trailing meniscus than at the leading meniscus, causing an increased flattening of the trailing end.
In Figure 4, the total drag acting on the bubble once it attains terminal velocity is shown for both the Newtonian and non-Newtonian cases. The drag coefficients are obtained by balancing the terminal-state total drag with the net buoyancy force, which yields . The total drag increases for increasing values of r0/R because of the increasing effects of occlusion and wall effects.34 As evident, the drag acting on the bubble for the non-Newtonian case is significantly higher than that acting on the Newtonian counterpart: the bubble moves more slowly in the non-Newtonian than in the Newtonian fluid. An important quantity of physiological significance is the wall shear stress and its variation at any given point on the tube wall during the steady-state motion of the bubble. The value of the wall shear stress for a fully developed channel flow of the non-Newtonian fluid in the absence of a bubble is denoted by τref. Again, for consistency, shear stress values for both the models are nondimensionalized with respect to τref.
In Figure 5, the time variation of shear stress acting on a given point on the tube wall is displayed for three different aspect ratios of the bubble for both the Newtonian and non-Newtonian models. t* refers to the nondimensional time (see Section 2 for nondimensional aspects). There is a reduced shear for non-Newtonian fluid (nearly half) compared with that of the Newtonian fluid. The reduction in wall shear is directly attributable to the decreased steady-state velocities of the bubble. The general features of the time variation of shear stress for both the Newtonian and non-Newtonian fluids are as follows. At initial times, the presence of the bubble is negligibly felt and the shear stress values are at base levels. Base value refers to the wall shear stress in a fully developed channel flow in the absence of a bubble. This is followed by a gradual increase in shear stress but in a direction opposite that of the base value. This is because as the bubble moves it entrains fluid to satisfy mass conservation. This entrainment is in a direction opposite that of the bulk flow. The entrained fluid is squeezed in the narrow gap between the bubble and the wall, giving rise to large wall shear stress acting in the opposite direction of the base value. After long times, the shear stress once again comes back to the base value, indicating the negligible effect of the bubble.
The time over which the wall shear stress varies from its base value will herein be referred to as “bubble residence time” because this is the period over which the presence of the bubble is felt. Figure 6 shows the bubble residence times for both the Newtonian and non-Newtonian cases. A smooth curve is generated by interpolating the values for various r0/R with known values at 0.5, 0.7, and 0.9. For Figure 8 and Figure 10, a similar procedure is carried out to obtain a smooth curve. The residence time increases with increasing aspect ratio. The values for the non-Newtonian case are almost twice that of the Newtonian case—that is, the effect of the bubble is felt longer for non-Newtonian fluid than for Newtonian fluid. Such elevated wall shear stress values coupled with increased residence time could result in an increased physiological response of the endothelial cells lining the vessel wall. For example, an increased shear stress may have significant effects on the modulation of cytosolic-free calcium (Ca2+) transients.38,39 For r0/R > 1, the steady-state velocity of the bubble is expected to asymptote to a constant value.19 Beyond r0/R > 1, the residence time will thus depend linearly on the length of the bubble. For r0/R > 1, assuming the shape of the bubble to be hemispherical capped cylinder, the length of the cylindrical bubble is proportional to (r0/R)3. Thus, the residence time is also expected to be proportional to (r0/R)3.
This case corresponds to a high shear rate condition (Umax/R = 180), where blood behaves almost like a Newtonian fluid of viscosity, μ∞. The size of the tube (diameter) is taken to be 1250 µm. The size of the plasma layer (ε) is assumed to be zero and the entire bulk fluid is modeled as a shear-thinning generalized power-law fluid. The nondimensional value of dp/dz = 0.85 and the same value is used for comparison with Newtonian fluid of viscosity, μ∞.
In Figure 7, the steady-state bubble shapes are plotted for aspect ratios of 0.5, 0.7, and 0.9. The vertical sidewalls of the bounding boxes denote the positions of the actual walls of the vertical tube. The bubble remains nearly spherical for all aspect ratios. Occlusion of the bubble in the tube is more pronounced for r0/R = 0.9.
In Figure 8, the steady-state total drag coefficient of the bubble is plotted against the aspect ratio of the bubble for the both the Newtonian and non-Newtonian cases. The drag coefficient steeply increases as the aspect ratio increases beyond 0.7 for both cases. The values of the drag coefficients are nearly identical in both cases, corroborating the fact that under high shear rates blood behavior is Newtonian.
In Figure 9, the time variation of shear stress at a given point on the tube wall is plotted similar to Figure 5. As expected, the values and the variations of shear stresses are identical for both cases for all aspect ratios. Complementary to Figure 9, the bubble residence time is plotted as a function of the aspect ratio in Figure 10. The residence times are identical for both the cases and a steep increase for aspect ratios beyond 0.7 is noteworthy. As stated earlier, this might have significant effect on the modulation of cytosolic-free calcium (Ca2+) transients for the endothelial cells that line the tube wall.
The flow features and the associated dynamics of a nearly occluding bubble in a shear-thinning, non-Newtonian blood were studied under two different flow conditions. These flow conditions corresponded to two widely diverse shear rates, Reynolds numbers, capillary numbers, Weber numbers, and Froude numbers. The rheology of blood was described using a generalized power-law model. A two-fluid model accounting for the cell-depleted or the plasma layer was used, where applicable. Results were then compared to a fully Newtonian model assuming blood viscosity to be Newtonian. The associated wall shear stresses and bubble residence times were also evaluated in all the cases. Whereas the results for the Newtonian and the non-Newtonian models agree well at high shear rates, they show complete disagreement at lower shear rates. In particular, the bubble residence times at lower shear rates were almost twice as high in the non-Newtonian model as in the Newtonian case. This finding may significantly affect the determination of the fate of endothelial cells that line the walls of blood vessels.
Thus, to properly describe the bubble dynamics for all shear rates prevalent in the blood vessels, a Newtonian blood viscosity model is inadequate and a non-Newtonian model is necessary.
We gratefully acknowledge NIH Grant R01 HL 67986-01 A1 and NASA Grant NNC05GA30G for supporting this study.
Conflicts of Interest
The authors declare no conflicts of interest.