Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
Ann N Y Acad Sci. Author manuscript; available in PMC 2010 April 1.
Published in final edited form as:
PMCID: PMC2790045

Bubble Motion through a Generalized Power-Law Fluid Flowing in a Vertical Tube


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.

Keywords: gas embolism, bubble, front tracking, shear-thinning fluid, power-law fluid, wall shear stress, bubble residence time

1. Introduction

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. 25 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.79 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. 1820). 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,2527 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. 3032, 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.

2. Mathematical Formulation

2.1 Fluid Mechanics: Formulation

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:

[nabla][center dot]u=0 and

ρ([partial differential]u[partial differential]t+[nabla]u[center dot]u) =[nabla]p+[nabla][center dot]μ([nabla]u+[nabla]Tu) +ρg+S(t)σκnfδ(xxf)ds.
Figure 1
Schematic of the problem.

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 δ(xxf) 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:

μl={ μp, Rεr<Rλ|γ˙|n1, 0r<Rε,



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 [gamma with dot above]· μ 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. [gamma with dot above] is related to the strain rate tensor (D) and is given as follows:

γ˙=2tr(D2), D=[nabla]u+[nabla]Tu2.

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.

2.1.1 Boundary Conditions

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:

[partial differential]u[partial differential]z(r,H,t)=0.

Symmetry boundary (r = 0): A reflective symmetry boundary condition is prescribed at r = 0,

[partial differential]uz[partial differential]r(0,H,t)=0,ur(0,H,t)=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.

2.1.2 Initial Condition

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.

2.2 Nondimensionalization Aspects

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 12ρlUmax2 max for pressure. With this scheme, the following nondimensional parameters that govern the motion are obtained:

Re=ρlUmax(2R)μ,Ca=μUmaxσWe=ρlUmax2(2R)σ=Re[center dot]Ca,Fr=Umax2gR.

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 [gamma with dot above] → ∞, μl → μ for the shear-thinning fluid. The steady state velocity of the bubble is represented by U in our subsequent discussions.

2.3 Numerical Methodology

We numerically solve Eqs. (1)(6) subject to boundary conditions given by Eqs. (7)(10). For completeness, we now briefly describe the numerical methodology.

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. 3537. 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 [gamma with dot above] 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 [gamma with dot above] (= [gamma with dot above] min) such that a stable convergent solution is obtained. [gamma with dot above] 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 [gamma with dot above]min = 0.2 resulted in stable converged solutions with negligible loss in accuracy.

3. Results and Discussion

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.

Figure 2
Re = 0.06: Comparison of inlet velocity profiles for non-Newtonian and Newtonian models for the same value of nondimensional dp/dz (=1100).

In surface tension–dominated regimes, which occur when values of Re [double less-than sign] 1, Ca [double less-than sign] 1 and for Re > 1, We [double less-than sign] 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.

3.1 Low Shear Rate: Re = 0.06

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.

Figure 3
Re = 0.06: Plot of steady-state bubble shapes for (i) r0/R = 0.5, (ii) r0/R = 0.7, and (iii) r0/R = 0.9.

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 (=CD(12ρlU2πr02)) with the net buoyancy force, which yields CD=8r0g3U2. 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.

Figure 4
Re = 0.06: Plot of steady-state bubble drag coefficient as a function of aspect ratio.

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.

Figure 5
Re = 0.06: Time variation of shear stress for a given point on the tube wall after the bubble moves with terminal velocity. |τref| denotes magnitude of wall shear stress for a fully developed flow of the non-Newtonian fluid in the absence of a ...

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.

Figure 6
Re = 0.06: Plot of bubble residence time versus aspect ratio. Residence time is denoted as the time over which there is an appreciable variation (~10%) of shear stress from the base value (absence of bubble) for any given point on the tube wall.
Figure 8
Re = 40: Plot of steady-state bubble drag coefficient as a function of aspect ratio.
Figure 10
Re = 40: Plot of bubble residence time versus aspect ratio. Residence time is denoted as the time over which there is an appreciable variation (~10%) of shear stress from the base value (absence of bubble) for any given point on the tube wall.

3.2 High Shear Rate: Re = 40

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.

Figure 7
Re = 40: Plot of steady-state bubble shapes for (i) r0/R = 0.5, (ii) r0/R = 0.7, and (iii) 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.

Figure 9
Re = 40: Time variation of shear stress for a given point on the tube wall after the bubble moves with terminal velocity. |τref| denotes magnitude of wall shear stress for a fully developed flow of the non-Newtonian fluid in the absence of a bubble. ...

4. Conclusions

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.


surface tension
curvature on the interface
Dirac delta function
consistency coefficient
power-law index
Wall shear stress
Plasma layer thickness
[gamma with dot above]
local strain rate
Δμ, Δn
generalized power-law model parameters
a, b, c, d
generalized power-law model parameters
arc length parameter
dynamic pressure
r, R
axial position
Heaviside function
velocity vector
u, U
velocity components
normal vector
position vector
gravity vector
strain rate tensor
Reynolds number
capillary number
Weber number
Froude number
gas phase
liquid phase
front point
infinite strain rate
residence time
time level


Conflicts of Interest

The authors declare no conflicts of interest.


1. Souders JE, et al. Spatial distribution of venous gas emboli in the lungs. J. Appl. Physiol. 1999;87:1937–1947. [PubMed]
2. Cavanagh DP, Eckmann DM. Interfacial dynamics of stationary gas bubbles in flows in inclined tubes. J. Fluid Mech. 1999;398:225–244.
3. Eckmann DM, Lomivorotov VN. Microvascular gas embolization clearance following per-fluorocarbon administration. J. Appl. Physiol. 2003;94:860–868. [PubMed]
4. Suzuki A, Eckmann DM. Embolism bubble adhesion force in excised perfused microvessels. Anesthesiology. 2003;99:400–408. [PubMed]
5. Zhang J, Eckmann DM, Ayyaswamy PS. A front tracking method for a deformable intravascular axisymmetric bubble with a soluble surfactant. J. Comp. Phys. 2006;214:366–396.
6. Popel AS, Johnson PC. Microcirculation and hemorheology. Annu. Rev. Fluid Mech. 2005;37:43–69. [PMC free article] [PubMed]
7. Cho YI, Kensey KR. Effects of the non-Newtonian viscosity of blood on flows in a diseased arterial vessel. 1. Steady flows. Biorheology. 1991;28:241–262. [PubMed]
8. Johnston BM, et al. Non-Newtonian blood flow in human right coronary arteries: steady state simulations. J. Biomech. 2004;37:709–720. [PubMed]
9. Ballyk PD, Steinman DA, Ethier CR. Simulation of non-Newtonian blood-flow in an end-to-side anastomosis. Biorheology. 1994;31:565–586. [PubMed]
10. Pries AR, Neuhaus D, Gaehtgens P. Blood-viscosity in tube flow—dependence on diameter and hematocrit. Am. J. Physiol. 1992;263:H1770–H1778. [PubMed]
11. Das B, Johnson PC, Popel AS. Computational fluid dynamic studies of leukocyte adhesion effects on non-Newtonian blood flow through microvessels. Biorheology. 2000;37:239–258. [PubMed]
12. Sharan M, Popel AS. A two-phase model for flow of blood in narrow tubes with increased effective viscosity near the wall. Biorheology. 2001;38:415–428. [PubMed]
13. Chhabra RP. Bubbles, Drops and Particles in Non-Newtonian Fluids. Boca Raton, FL: CRC Press; 1993.
14. Rodrigue D, De Kee D, Fong CFCM. The slow motion of a single gas bubble in a non-Newtonian fluid containing surfactants. J. Non-Newtonian Fluid Mech. 1999;86:211–227.
15. Rodrigue D. A simple correlation for gas bubbles rising in power-law fluids. Can. J. Chem. Eng. 2002;80:289–292.
16. Tsukada T, et al. Theoretical and experimental studies of the deformation of bubbles moving in quiescent Newtonian and non-Newtonian liquids. J. Chem. Eng. Jpn. 1990;23:192–198.
17. Ohta M, et al. Dynamic processes in a deformed drop rising through shear-thinning fluids. J. Non-Newtonian Fluid Mech. 2005;132:100–107.
18. Hyman WA, Skalak R. Viscous flow of a suspension of liquid drops in a cylindrical tube. Appl. Sci. Res. 1972;26:27–52.
19. Martinez MJ, Udell KS. Axisymmetric creeping motion of drops through circular tubes. J. Fluid Mech. 1990;210:565–591.
20. Tsai TM, Miksis MJ. Dynamics of a drop in a constricted capillary tube. J. Fluid Mech. 1994;274:197–217.
21. Tryggvason G, et al. A front-tracking method for the computations of multiphase flow. J. Comp. Phys. 2001;169:708–759.
22. Osher S, Fedkiw RP. Level set methods: an overview and some recent results. J. Comput. Phys. 2001;169:463–502.
23. Jacqmin D. Calculation of two-phase Navier-Stokes flows using phase-field modeling. J. Comput. Phys. 1999;155:96–127.
24. Yue PT, et al. Phase-field simulations of interfacial dynamics in viscoelastic fluids using finite elements with adaptive meshing. J. Comput. Phys. 2006;219:47–67.
25. Hirt CW, Nichols BD. Volume of fluid (VOF) methods for the dynamics of free boundaries. J. Comput. Phys. 1981;39:201–225.
26. Scardovelli R, Zaleski S. Direct numerical simulation of free-surface and interfacial flow. Annu. Rev. Fluid Mech. 1999;31:567–603.
27. James AJ, Lowengrub J. A surfactant-conserving volume-of-fluid method for interfacial flows with insoluble surfactant. J. Comput. Phys. 2004;201:685–722.
28. Sussman M, Puckett EG. A coupled level set and volume-of-fluid method for computing 3D and axisymmetric incompressible two-phase flows. J. Comput. Phys. 2000;162:301–337.
29. LeVeque RJ, Li Z. Immersed interface methods for Stokes flow with elastic boundaries or surface tension. SIAM J. Sci. Statist. Comput. 1997;18:709–735.
30. Fedkiw RP, et al. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method) J. Comput. Phys. 1999;152:457–492.
31. Udaykumar HS, et al. An Eulerian method for computation of multimaterial impact with ENO shock-capturing and sharp interfaces. J. Comput. Phys. 2003;186:136–177.
32. Quan SP, Schmidt DP. A moving mesh interface tracking method for 3D incompressible two-phase flows. J. Comput. Phys. 2007;221:761–780.
33. Shin S, et al. Accurate representation of surface tension using the level contour reconstruction method. J. Comput. Phys. 2005;203:493–516.
34. Mukundakrishnan K, et al. Numerical study of wall effects on buoyant gas-bubble rise in a liquid-filled finite cylinder. Phys. Rev. E Stat. Nonlin. Soft Matter Phys. 2007;76 036308. [PMC free article] [PubMed]
35. Almgren AS, et al. A conservative adaptive projection method for the variable density incompressible Navier-Stokes equations. J. Comp. Phys. 1998;142:1–46.
36. Bell JB, Colella P, Glaz HM. A 2nd-order projection method for the incompressible Navier Stokes equations. J. Comput. Phys. 1989;85:257–283.
37. Puckett EG, et al. A high-order projection method for tracking fluid interfaces in variable density incompressible flows. J. Comput. Phys. 1997;130:269–282.
38. Wiesner TF, Berk BC, Nerem RM. A mathematical model of cytosolic calcium dynamics in human umbilical vein endothelial cells. Am. J. Physiol. Cell Physiol. 1996;39:C1556–C1569. [PubMed]
39. Wiesner TF, Berk BC, Nerem RM. A mathematical model of the cytosolic-free calcium response in endothelial cells to fluid shear stress. Proc. Natl. Acad. Sci. USA. 1997;94:3726–3731. [PubMed]