Flow Capillary Number Research Articles

On the extension of Bretherton theory for thin liquid films formed around elongated bubbles

The formation of thin liquid films around an elongated bubble moving in a capillary is pertinent to many applications. However, development of a theoretical model for the film thickness has been a challenge for several decades. The prominent theory characterizing the film thickness was developed by Bretherton. This theory relates the liquid film thickness in axisymmetric capillaries to the flow capillary number (Ca) for values of up to 0.003. Modified forms of the Bretherton theory have been presented for Ca values as high as 2. However, the validity of these models has not been rigorously examined. In addition, the validity of the Bretherton model itself in non-axisymmetric cross section capillaries remains uncharted. The objective of this paper is to determine whether the Bretherton relation can be extended to a broader range of Ca values and non-axisymmetric cross section channels, wherein the film thickness is not uniform along the channel circumference. A series of experiments are conducted on a set of fluids and different channel sizes and cross-sectional geometries to produce a wide range of viscous, surface tension, and inertial forces. The results show that when inertial forces are significant, and modified Bretherton models fail to predict the film thickness at Ca well below 2. The experimental results also show that depending on fluid properties, another key requirement of Bretherton solution of Landau–Levich equation may not be met. In rectangular cross section channels, the difference in surface tension forces along the longer and shorter channel axes, results in great deviation from the existing lubrication-based theories. This deviation greatly expands with the increase in the channel aspect ratio.

On methods to reduce spurious currents within VOF solver frameworks. Part 1: a review of the static bubble/droplet

Abstract In two-phase flows in which the Capillary number is low, errors in the computation of the surface tension force at the interface cause Front-Capturing methods such as Volume of Fluid (VOF) and Level-Set (LS) to develop interfacial spurious currents. To better solve low Capillary number flows, special treatment is required to reduce such spurious currents. Smoothing the phase indicator field to more accurately compute the curvature or adding interfacial artificial viscosity are techniques that can treat this problem. This study explores OpenFOAM, Fluent and StarCCM+ VOF solvers for the classical case of a static bubble/droplet immersed in a continuous aqueous phase, with the focus on the ability of these solvers to adequately reduce spurious currents. The results are expected to be helpful for practicing chemical engineers who use multiphase CFD solvers in their work.

Thin liquid films formation and evaporation mechanisms around elongated bubbles in rectangular cross-section microchannels

Thin film evaporation is the dominant mechanism of heat transfer in the microchannel flow boiling process. Hence, its understanding is paramount to development of predictive models and more efficient microchannel heat exchangers. Although microchannels mostly have a rectangular cross-section, the existing models are developed for circular cross-section microchannels. The non-uniformity of the liquid film thickness in microchannels with a rectangular cross-section has made prediction of the onset of dryout and heat transfer coefficient greatly challenging. Here, a test device capable of measuring the wall heat flux is utilized to fully characterize the thin film formation and evaporation process. The effects of flow capillary number (Ca) and channel aspect ratio on the liquid film thickness are determined in microchannels with 300 μm width and 300, 150, and 75 μm heights, representing aspect ratios of 1, 2, and 4, respectively. It is shown that the only existing mechanistic model for the elongated bubble regime in microchannels fails to predict the experimental heat transfer coefficient (HTC) both quantitively, by a 50–120% margin, and qualitatively, due to a fundamental difference in evaporation and dryout of a variable versus a uniform film assumed in the model. The results presented here pave the way for development of a more accurate mechanistic model for the thin film evaporation heat transfer mechanism in microchannels with a rectangular cross-section.

Computations of coating windows for reverse roll coating of liquid films

This paper focuses on the flow stability of reverse roll coating flows employed in practical applications to coat thin liquid films on substrates of different materials. Other coating methods are available and various computer programs are purported to predict (for each corresponding method) the set of parameters that permit stable coating conditions. The application of these programs yields so-called coating windows for each particular coating setup. Knowing these coating windows permits coating methods to be applied with preknowledge of the sets of parameters for stable coatings. For reverse roll coating, computer codes for predicting coating windows are scarce or even nonexistent, which was the reason for the present work on the analytical considerations described in this paper. The results obtained formed the basis for a reverse roll coating computer program and its applications. The program is described, and computational results are provided. These are compared with experimental data available in the literature. Good agreement was obtained for high capillary number flows, but there were discrepancies between experimental and computational results for lower capillary numbers. The experimental results for low capillary number flows also appeared less reliable, which triggered the present experimental work. Some experimental results are given for low capillary numbers and these results are also compared with corresponding numerical data. It is suggested that a test section should be built consisting of two rollers and one slot die to study experimentally the stability regions of reverse roll coating flows in a rigorous manner.

Dynamic pore-network modeling of air-water flow through thin porous layers

Thin porous layers, that have large aspect ratios, are seen in many applications such as hydrogen fuel cells and hygiene products, in which air-water immiscible flow is of great interest. Direct numerical simulations based on Navier-Stokes equation are computationally expensive, and even prohibitive for low capillary number flow such as water flooding in low-temperature polymer electrolyte fuel cells. Alternatively, the pore-network modeling needs much less computational resources, while still retaining essentials of the pore-structure information. In this work, a dynamic pore-network model of air-water flow with phase change has been developed. We focus on drainage processes through thin porous layers, in which liquid water is the nonwetting phase. Three test cases are conducted, namely, air-water flow through a thin porous layer, air-water flow through a bilayer of fine and coarse thin porous layers, and water flooding in the gas diffusion layer of a polymer electrolyte fuel cell with phase change between water and its vapor. Using these test cases, we aim to demonstrate the application of dynamic pore-network modeling in thin porous media studies. In particular, we discuss the challenge of modeling thin porous media at the average scale, and highlight the role of phase change in removing liquid water from the cathode gas diffusion layer.

A CutFEM method for two-phase flow problems

In this article, we present a cut finite element method for two-phase Navier–Stokes flows. The main feature of the method is the formulation of a unified continuous interior penalty stabilisation approach for the stabilisation of the advection and the pressure–velocity coupling, on the one hand, and the stabilisation of the cut region on the other hand. The accuracy of the algorithm is enhanced by the development of extended fictitious domains to guarantee a well defined velocity from previous time steps in the current domains. Finally, the robustness of the moving-interface algorithm is further improved by the introduction of a curvature smoothing technique that reduces spurious velocities. The algorithm is shown to perform remarkably well for low capillary number flows, and is a first step towards flexible and robust CutFEM algorithms for the simulation of flows in microfluidic devices.

Dewetting of Thin Liquid Films Surrounding Air Bubbles in Microchannels.

As an air bubble translates in a microchannel, a thin film of liquid is formed on the bounding walls. In a microchannel with a rectangular cross-section, the liquid in the film leaks toward the low-pressure corners of the geometry, which leads to the appearance of local minima in the film thickness in the cross-sectional plane. In such a configuration, theory suggests that the minimum film thickness scales with Ca and Ca4/3 depending on the distance from the nose of the bubble, where Ca = μUb/γ is the flow capillary number based on the bubble velocity Ub, liquid viscosity μ, and surface tension γ, and Ca ≪ 1. We show that the film of a partially wetting liquid dewets on the channel wall at the sites of the local minima in the film thickness as it acquires thicknesses around and below 100 nm. Our experiments show that the distance Lw between the nose of the bubble and the initial dewetting location is a function of Ca and surface wettability. For channels of different wettability, Lw always scales proportional to Caα, where 1.7 < α < 2 for the range of 10-5 < Ca < 10-2. Moreover, Lw increases up to 10 times by enhancing the wettability of the surface at a given Ca. Our present measurements of Lw provide a design constraint on the lengths of bubbles to maintain a liquid wet channel without dry patches on the wall.

Effect of wettability on two-phase quasi-static displacement: Validation of two pore scale modeling approaches

Understanding of pore-scale physics for multiphase flow in porous media is essential for accurate description of various flow phenomena. In particular, capillarity and wettability strongly influence capillary pressure-saturation and relative permeability relationships. Wettability is quantified by the contact angle of the fluid-fluid interface at the pore walls. In this work we focus on the non-trivial interface equilibria in presence of non-neutral wetting and complex geometries. We quantify the accuracy of a volume-of-fluid (VOF) formulation, implemented in a popular open-source computational fluid dynamics code, compared with a new formulation of a level set (LS) method, specifically developed for quasi-static capillarity-dominated displacement. The methods are tested in rhomboidal packings of spheres for a range of contact angles and for different rhomboidal configurations and the accuracy is evaluated against the semi-analytical solutions obtained by Mason and Morrow (1994). While the VOF method is implemented in a general purpose code that solves the full Navier-Stokes (NS) dynamics in a finite volume formulation, with additional terms to model surface tension, the LS method is optimized for the quasi-static case and, therefore, less computationally expensive. To overcome the shortcomings of the finite volume NS-VOF system for low capillary number flows, and its computational cost, we introduce an overdamped dynamics and a local time stepping to speed up the convergence to the steady state, for every given imposed pressure gradient (and therefore saturation condition). Despite these modifications, the methods fundamentally differ in the way they capture the interface, as well as in the number of equations solved and in the way the mean curvature (or equivalently capillary pressure) is computed. This study is intended to provide a rigorous validation study and gives important indications on the errors committed by these methods in solving more complex geometry and dynamics, where usually many sources of errors are interplaying.

Microbubble moving in blood flow in microchannels: effect on the cell-free layer and cell local concentration

Gas embolisms can hinder blood flow and lead to occlusion of the vessels and ischemia. Bubbles in microvessels circulate as tubular bubbles (Taylor bubbles) and can be trapped, blocking the normal flow of blood. To understand how Taylor bubbles flow in microcirculation, in particular, how bubbles disturb the blood flow at the scale of blood cells, experiments were performed in microchannels at a low Capillary number. Bubbles moving with a stream of in vitro blood were filmed with the help of a high-speed camera. Cell-free layers (CFLs) were observed downstream of the bubble, near the microchannel walls and along the centerline, and their thicknesses were quantified. Upstream to the bubble, the cell concentration is higher and CFLs are less clear. While just upstream of the bubble the maximum RBC concentration happens at positions closest to the wall, downstream the maximum is in an intermediate region between the centerline and the wall. Bubbles within microchannels promote complex spatio-temporal variations of the CFL thickness along the microchannel with significant relevance for local rheology and transport processes. The phenomenon is explained by the flow pattern characteristic of low Capillary number flows. Spatio-temporal variations of blood rheology may have an important role in bubble trapping and dislodging.

Dynamics of isolated confined air bubbles in liquid flows through circular microchannels: an experimental and numerical study

Experimental and numerical studies are performed to characterise the dynamics of isolated confined air bubbles in laminar fully developed liquid flows within channels of diameters d = 0.5 mm and d = 1 mm. Water and glycerol are used as the continuous liquid phase, and therefore, a large range of flow capillary numbers 10−4 < Ca < 10−1 and Reynolds numbers 10−3 < Re < 103 are covered. An extensive investigation is performed on the effect of bubble size and flow capillary number on different flow parameters, such as the shape and velocity of bubbles, thickness of the liquid film formed between the bubbles and the channel wall, and the development lengths in front and at the back of the bubbles. The micro-particle shadow velocimetry technique (μPSV) is employed in the experimental measurements allowing simultaneous quantification of important flow parameters using a single sequence of high-speed greyscale images recorded at each test condition. Bubble volume and flow rate of the continuous liquid phase are precisely determined in the post-processing stage using the μPSV images. These parameters are then used as initial and boundary conditions to set up CFD simulations reproducing the corresponding two-phase flow. Simulations based on the volume of fluid technique with the aim of capturing the interface dynamics are performed with both ANSYS Fluent v. 14.5, here augmented by implementing self-defined functions to improve the accuracy of the surface tension force estimation, and ESI OpenFOAM v. 2.1.1. The present approach not only results in valuable findings on the underlying physics involved in the problem of interest but also allows us to directly compare and validate results that are currently obtained by the experimental and computational methods. It is believed that similar methodology can be employed to rigorously investigate more complex two-phase flow regimes in micro-geometries.

Characterization and prediction of CO2 saturation resulting from modeling buoyant fluid migration in 2D heterogeneous geologic fabrics

Understanding the pore saturation of CO2 injected into heterogeneous rocks for permanent storage is a key challenge for constraining storage efficiency (capacity), plume migration extent, and seismic imaging. Much work on fluid flow and saturation has focused on traditional pressure-gradient driven flow represented by viscous multi-phase Darcy flow. Here we investigate CO2 saturations resulting from buoyancy-driven capillary flow, where buoyancy forces dominate viscous forces. These low capillary number flow regimes (Ca<10−4) exist some distance from injection wells, potentially representing the majority of the storage domain during and after injection. We simulate CO2-brine buoyant displacement patterns using invasion percolation (IP) methods in various decimeter-scale, highly resolved heterogeneous 2D models. Different sedimentologic fabrics result in dramatic variability of total CO2 saturation. We seek a generalized predictive model that allows CO2 saturation resulting from capillary flow to be estimated reasonably for 2D domains from fundamental geologic and fluid properties. We focus on three metrics of pore threshold pressure in clastic porous materials: the central value (median; mean), the range (standard deviation), and the ratio of the horizontal to vertical correlation length. These are formalized in the description of 54 clastic facies and three fabrics. An ‘IP Saturation Predictor’ model is derived by linear interpolation of simulation results from synthetic stochastic models representing a wide range of facies and fabrics. To validate model results, we compare the predictive model with flow simulations in model domains extracted from highly resolved natural geologic sedimentary samples at similar resolution.

THERMAL TRANSPORT BEHAVIOR OF A LIQUID PLUG MOVING INSIDE A DRY CAPILLARY TUBE

Computational simulation of flow and heat -transfer characteristics of a single isolated liquid plug (of glycerin and water, respectively) flowing in a round capillary tube (internal diameter, D = 2.0 mm), so as to understand its local thermo-hy- drodynamic transport behavior, is reported. Both, constant heat flux and constant-tem - perature boundary conditions have been considered at the tube wall. When an isolated liquid plug moves in a capillary tube, there is a difference in the advancing and receding dynamic contact angles of the two menisci, respectively. This has been considered in the simulations. The linearized simplification of Hoffman-Tanner's law is used to model the variation in the two respective apparent dynamic contact angles with the velocity of the liquid plug (i.e., Ca = µU/σ). Simulations are carried out for a range of Capillary numbers and lengths of liquid plugs (L/D ratio). It has been found that variation in dy- namic contact angle leads to an enhanced local and average heat-transfer coefficient in the moving liquid plug; the local fluid circulation being affected by menisci deformation. In addition, as the length of the liquid plug is increased, the heat-transfer coefficient decreases and finally shows the asymptotic transport behavior of Poiseuille flow. Other than the flow Capillary number, the fluid Prandtl number and the L/D ratio of the liquid plug, strongly affect the local thermo hydrodynamics.

Shear-induced particle migration and margination in a cellular suspension

We simulate the cross-flow migration of rigid particles such as platelets in a red blood cell (RBC) suspension using the Stokes flow boundary integral equation method. Two types of flow environments are investigated: a suspension undergoing a bulk shear motion and a suspension flowing in a microchannel or duct. In a cellular suspension undergoing bulk shear deformation, the cross-flow migration of particles is diffusional. The velocity fluctuations in the suspension, which are the root cause of particle migration, are analyzed in detail, including their magnitude, the autocorrelation of Lagrangian tracer points and particles, and the associated integral time scales. The orientation and morphology of red blood cells vary with the shear rate, and these in turn cause the dimensionless particle diffusivity to vary non-monotonically with the flow capillary number. By simulating RBCs and platelets flowing in a microchannel of 34 μm height, we demonstrate that the velocity fluctuations in the core cellular flow region cause the platelets to migrate diffusively in the wall normal direction. A mean lateral velocity of particles, which is most significant near the edge of the cell-free layer, further expels them toward the wall, leading to their excess concentration in the cell-free layer. The calculated shear-induced particle diffusivity in the cell-laden region is in qualitative agreement with the experimental measurements of micron-sized beads in a cylindrical tube of a comparable diameter. In a smaller duct of 10 × 15 μm cross section, the volume exclusion becomes the dominant mechanism for particle margination, which occurs at a much shorter time scale than the migration in the bigger channel.

Numerical Simulation of the Formation of a Jet and Droplets from a Capillary in a Co-Flowing Ambient Fluid

Formation of a liquid jet and droplets from a micro-capillary in another immiscible co-flowing ambient fluid is simulated numerically with a front-tracking/finite difference method. The transition from dripping to jetting is observed and the size of the resulting droplets is predicted. The predicted size in jetting mode is in close agreement with that from linear stability theory for a moving liquid jet in a co-flowing ambient fluid with the same speed as the jet. The transition between dripping and jetting are mapped on the diagram of the jet Weber number versus the ambient flow Capillary number. It is found that the viscosity ratio is an important parameter to distinguish the mode.

Dynamics of viscoelastic liquid filaments: Low capillary number flows

Many applications of viscoelastic free surface flows requiring formation of drops from small nozzles, e.g., ink-jet printing, micro-arraying, and atomization, involve predominantly extensional deformations of liquid filaments. The capillary number, which represents the ratio of viscous to surface tension forces, is small in such processes when drops of water-like liquids are formed. The dynamics of extensional deformations of viscoelastic liquids that are weakly strain hardening, i.e., liquids for which the growth in the extensional viscosity is small and bounded, are here modeled by the Giesekus, FENE-P, and FENE-CR constitutive relations and studied at low capillary numbers using full 2D numerical computations. A new computational algorithm using the general conformation tensor based constitutive equation [M. Pasquali, L.E. Scriven, Theoretical modeling of microstructured liquids: a simple thermodynamic approach, J. Non-Newtonian Fluid Mech. 120 (2004) 101–135] to compute the time dependent viscoelastic free surface flows is presented. DEVSS-TG/SUPG mixed finite element method [M. Pasquali, L.E. Scriven, Free surface flows of polymer solutions with models based on conformation tensor, J. Non-Newtonian Fluid Mech. 108 (2002) 363–409] is used for the spatial discretization and a fully implicit second-order predictor–corrector scheme is used for the time integration. Inertia, capillarity, and viscoelasticity are incorporated in the computations and the free surface shapes are computed along with all the other field variables in a fully coupled way. Among the three models, Giesekus filaments show the most drastic thinning in the low capillary number regime. The dependence of the transient Trouton ratio on the capillary number in the Giesekus model is demonstrated. The elastic unloading near the end plates is investigated using both kinematic [M. Yao, G.H. McKinley, B. Debbaut, Extensional deformation, stress relaxation and necking failure of viscoelastic filaments, J. Non-Newtonian Fluid Mech. 79 (1998) 469–501] and energy analyses. The magnitude of elastic unloading, which increases with growing elasticity, is shown to be the largest for Giesekus filaments, thereby suggesting that necking and elastic unloading are related.

Capillary Filling Flows inside Patterned‐Surface Microchannels

AbstractCapillary flows inside microchannels with patterned‐surfaces are investigated theoretically and numerically. The surface energy method is used to derive an equivalent contact angle (ECA) model for small capillary number flows. The SIMPLE algorithm using a volume of fluid (VOF) method is adopted to investigate the flows in those microchannels. The flow characteristics such as the liquid front shapes and the evolution of the liquid lengths are obtained. The numerical results reveal that capillary flows in a patterned‐surface microchannel still follow the traditional capillary theories. The ECA model is confirmed by the numerical results. It indicates that the capillary flows inside the patterned‐surface microchannels can be estimated by means of the homogeneous‐surface microchannels with the equivalent contact angle. The ECA model provides a good criterion for the total wettability of a patterned‐surface microchannel, as well.

Gas displacing liquids from non-circular tubes: high capillary number flow of a shear-thinning liquid

Transient, three-dimensional finite element analysis has been used to investigate the displacement of a shear thinning liquid from prismatic channels of square, rectangular and trapezoidal cross sections. Inertia, gravity and surface tension effects are neglected and the results therefore apply in the limits of low Reynolds and high capillary numbers. The analysis is carried out in a fixed frame of reference and gas penetration is modelled as the bubble moves down the tube, which is long relative to its transverse dimension. Results are provided for the thickness of the layers left on the channel walls under developed conditions, and the fraction of the cross section occupied by liquid, as a function of the channel cross-sectional geometry and the degree of shear thinning, modelled using the power law. Interface contours on the channel cross sections are displayed. It is found for the Newtonian liquid that the fingering instability arises in the rectangular channel when the aspect ratio reaches about five. Shear thinning delays the onset of the instability to higher aspect ratios. The results are systematized, and insights gained into the influence of channel geometry and shear thinning, by noting a qualitative, inverse relationship between the deposited layer thickness and the shear rate at the wall in the flow ahead of the bubble.

Gas displacing liquids from tubes: high capillary number flow of a power law liquid including inertia effects

Computer simulation has been used as a virtual experimental tool to investigate the displacement of a shear thinning Power Law liquid from a cylindrical tube by a gas, in the limit of high capillary number and in the absence of gravity effects. Two scenarios have been considered. In the first, gas enters at a steady rate, and the gas penetration velocity and residual wall layer thickness attain steady values. In the second, a constant gas pressure is applied at the inlet, and the gas penetration rate accelerates as the column of liquid ahead of it becomes shorter. The first set of experiments confirm that the developed wall layer thickness falls with increased degrees of shear thinning and with increased Reynolds Number, and quantifies the latter effect for the first time. The relationship is summarized by a correlation formula for dimensionless layer thickness as a function of Power Law index, n, and an appropriately defined Reynolds group in the range 0.1⩽ n⩽1.0,0.001⩽ Re⩽100. The flow pattern ahead of the gas bubble throughout the range of these experiments was always of the ‘by-pass’ type, consistent with a generalized criterion for the transition between by-pass and re-circulating flow which is derived for a Power Law liquid. In the second set of experiments, where a constant gas inlet pressure is applied, giving accelerating gas penetration, a comparison of layer thickness values at various axial positions, with those obtained at corresponding Reynolds number in steady flow, showed close agreement, though a small discrepancy for the highest Reynolds numbers could indicate some influence of inertia in the accelerating liquid column. At higher Reynolds number, in both steady and accelerating flow, the gas bubble near to the inlet shows a concave region on the axis, with re-circulation in the liquid ahead of it. As the bubble moves down the tube, the radius of this concavity decreases and a steady convex profile is eventually attained, with reversion of the flow to by-pass type. We show that the origin of this is inertial. The results have applications in a number of technologies, including gas-assisted injection moulding of plastics and certain gas liquid reactors.

Contact-line dynamics of a diffuse fluid interface

An investigation is made into the moving contact line dynamics of a Cahn–Hilliard–van der Waals (CHW) diffuse mean-field interface. The interface separates two incompressible viscous fluids and can evolve either through convection or through diffusion driven by chemical potential gradients. The purpose of this paper is to show how the CHW moving contact line compares to the classical sharp interface contact line. It therefore discusses the asymptotics of the CHW contact line velocity and chemical potential fields as the interface thickness ε and the mobility κ both go to zero. The CHW and classical velocity fields have the same outer behaviour but can have very different inner behaviours and physics. In the CHW model, wall–liquid bonds are broken by chemical potential gradients instead of by shear and change of material at the wall is accomplished by diffusion rather than convection. The result is, mathematically at least, that the CHW moving contact line can exist even with no-slip conditions for the velocity. The relevance and realism or lack thereof of this is considered through the course of the paper.The two contacting fluids are assumed to be Newtonian and, to a first approximation, to obey the no-slip condition. The analysis is linear. For simplicity most of the analysis and results are for a 90° contact angle and for the fluids having equal dynamic viscosity μ and mobility κ. There are two regions of flow. To leading order the outer-region velocity field is the same as for sharp interfaces (flow field independent of r) while the chemical potential behaves like r−ξ, ξ = π/2/max{θeq, π − θeq}, θeq being the equilibrium contact angle. An exception to this occurs for θeq = 90°, when the chemical potential behaves like ln r/r. The diffusive and viscous contact line singularities implied by these outer solutions are resolved in the inner region through chemical diffusion. The length scale of the inner region is about 10√μκ – typically about 0.5–5 nm. Diffusive fluxes in this region are O(1). These counterbalance the effects of the velocity, which, because of the assumed no-slip boundary condition, fluxes material through the interface in a narrow boundary layer next to the wall.The asymptotic analysis is supplemented by both linearized and nonlinear finite difference calculations. These are made at two scales, experimental and nanoscale. The first set is done to show CHW interface behaviour and to test the qualitative applicability of the CHW model and its asymptotic theory to practical computations of experimental scale, nonlinear, low capillary number flows. The nanoscale calculations are carried out with realistic interface thicknesses and diffusivities and with various assumed levels of shear-induced slip. These are discussed in an attempt to evaluate the physical relevance of the CHW diffusive model. The various asymptotic and numerical results together indicate a potential usefullness for the CHW model for calculating and modelling wetting and dewetting flows.

Rheology of bubble-bearing magmas

The rheology of bubble-bearing suspensions is investigated through a series of three-dimensional boundary integral calculations in which the effects of bubble deformation, volume fraction, and shear rate are considered. The behaviour of bubbles in viscous flows is characterized by the capillary number, Ca, the ratio of viscous shear stresses that promote deformation to surface tension stresses that resist bubble deformation. Estimates of Ca in natural lava flows are highly variable, reflecting variations in shear rate and melt viscosity. In the low capillary number limit (e.g., in carbonatite flows) bubbles remain spherical and may contribute greater shear stress to the suspension than in high capillary number flows, in which bubble deformation is significant. At higher Ca, deformed bubbles become aligned in the direction of flow, and as a result, contribute less shear stress to the suspension. Calculations indicate that the effective shear viscosity of bubbly suspensions, at least for Ca<0.5, is a weakly increasing function of volume fraction and that suspensions of bubbles are shear thinning. Field observations and qualitative arguments, however, suggest that for sufficiently large Ca ( Ca greater than about 1) the effective shear viscosity may be less than that of the suspending liquid. Bubbles reach their quasi-steady deformed shapes after strains of order one; for shorter times, the continuous deformation of the bubbles results in continual changes of rheological properties. In particular, for small strains, the effective shear viscosity of the suspension may be less than that of the liquid phase, even for small Ca. Results of this study may help explain previous experimental, theoretical, and field based observations regarding the effects of bubbles on flow rheology.