Rheology of Dense Bubble Suspensions

The rheological behavior of rapidly sheared bubble suspensions is examined through numerical simulations and kinetic theory. The limiting case of spherical bubbles at large Reynolds number Re and small Weber number We is examined in detail. Here, Re (cid:53) (cid:114)(cid:103) a 2 / (cid:109) and We (cid:53) (cid:114)(cid:103) 2 a 3 / s , a being the bubble radius, (cid:103) the imposed shear, s the interfacial tension, and (cid:109) and (cid:114) , respectively, the viscosity and density of the liquid. The bubbles are assumed to undergo elastic bounces when they come into contact; coalescence can be prevented in practice by addition of salt or surface-active impurities. The numerical simulations account for the interactions among bubbles which are assumed to be dominated by the potential ﬂow of the liquid caused by the motion of the bubbles and the shear-induced collision of the bubbles. A kinetic theory based on Grad’s moment method is used to predict the distribution function for the bubble velocities and the stress in the suspension. The hydrodynamic interactions are incorporated in this theory only through their inﬂuence on the virtual mass and viscous dissipation in the suspension. It is shown that this theory provides reasonable predictions for the bubble-phase pressure and viscosity determined from simulations including the detailed potential ﬂow interactions. A striking result of this study is that the variance of the bubble velocity can become large compared with ( (cid:103) a ) 2 in the limit of large Reynolds number. This implies that the disperse-phase pressure and viscosity associated with the ﬂuctuating motion of the bubbles is quite signiﬁcant. To determine whether this prediction is reasonable even in the presence of nonlinear drag forces induced by bubble deformation, we perform simulations in which the bubbles are subject to an empirical drag law and show that the bubble velocity variance can be as large as


I. INTRODUCTION
The classic experiments of Bagnold 1 demonstrated that rapidly shearing a suspension of particles could induce both tangential and normal stresses that are much larger than those in the pure fluid. Bagnold's experiments were performed with neutrally buoyant solid particles suspended in a liquid, and the large stresses were observed when the Reynolds number and Stokes number of the particles were both large, indicating that both fluid and particle inertia were important. Subsequent experiments have shown that a similar behavior can be seen when dry granular materials are rapidly sheared. 2 A theoretical description of these effects has been developed for granular materials using a modification of the kinetic theory of dense gases that takes account of the inelasticity of the interparticle collisions. 3 In this theory, the tangential and normal stresses can be understood in terms of a particulate-phase effective viscosity and pressure. Sangani et al. 4 have extended this analysis to include the effects of a low Reynolds number flow of the interstitial gas. To date, however, there is no comparable theory for systems in which the inertia of the interstitial fluid plays an important role as it would when the continuous phase is liquid.
Theoretical analysis and numerical simulations in the presence of fluid inertia are generally quite difficult. How-ever, in the special case of a suspension of spherical bubbles at high Reynolds number, the fluid flow can be approximated as a potential flow. The equation for the velocity potential is then the linear, Laplace equation and detailed numerical simulations and kinetic theory can be developed. This approximation is valid if the bubbles' Reynolds number Reϭ␥a 2 / is large and their Weber number Weϭ␥ 2 a 3 /s is small. Here, s is the interfacial tension, and are the liquid viscosity and density, ␥ is the shear rate, and a is the bubble radius. These conditions can be achieved for a narrow range of bubble radii ͓with aϭO͑0.5 mm͔͒ in water. Throughout most of this paper, we will assume potential flow in the continuous phase and spherical bubbles because of the great theoretical simplifications that these approximations afford. However, in Section VI, we will consider the implications of the nonlinear drag force law arising at finite We for the qualitative predictions of the theory. In this paper, we will study rapid shear flows of bubble suspensions using numerical simulations that include the effects of the potential-flow interactions among the bubbles. We will also develop a kinetic theory based on Grad's moment method. 5 It will be seen that the hydrodynamic interactions can be incorporated in such a theory by a simple adjustment of the virtual mass and drag coefficient. As a result, the kinetic theory for a sheared bubbly liquid will be shown to be quite similar to that derived previously for gas-solid suspensions. 4 Therefore, we shall see that the variance of the bubble velocity becomes quite large leading to a large disperse-phase pressure and viscosity analogous to that predicted and measured previously in granular flows.
We have previously presented limited simulation results for very dilute bubble suspensions in a proceedings paper. 6 There, we showed that dilute bubble suspensions can exhibit multiple steady states at the same shear rate. These consisted of a quenched state, in which most of the bubbles followed the motion of the liquid, and an ignited state, in which the variance of the bubble velocity was large compared with (␥a) 2 . A kinetic theory similar to that developed for dilute gas-solid suspensions by Tsao and Koch 7 was able to predict this multiplicity of steady states.
In this paper, we will consider higher bubble volume fractions that are typically characteristic of flows in which the disperse phase has a large effect on the suspension rheology. One of the motivations for studying sheared bubble suspensions is to assess the importance of mean velocity gradients in describing flows of bubble suspensions through closed conduits such as vertical pipes. The potential flow interactions among bubbles rising in an otherwise quiescent liquid tend to create bubble clusters. 8,9 This cluster formation is particularly severe when the magnitude of bubble-phase velocity variance is small. One role of the shear-induced bubble pressure investigated here may be to stabilize the homogeneous state of the suspension against this clustering instability. Of course, since we consider here only the case of mean shear in the absence of buoyancy forces producing a mean relative motion, the present study cannot prove this conjecture. The problem of determining a quantitative criterion for the stability of bubbly liquids in the presence of both mean relative motion and mean shear, however, is considerably more involved and is therefore left to a future investigation. Finally, we note that, in addition to its significance to the above problem, the shear-induced bubble pressure will also tend to prevent variations in the bubble volume fraction driven by other types of forces, such as lift and centrifugal forces, that may arise in pipe flows, vortical flows, and flows through horizontal channels. Although the present study is restricted to linear shear flows, the results obtained here will provide insight into these more complex flows.
In Section II, we review the basic equations governing the trajectories of spherical bubbles in a low viscosity liquid containing a sufficient amount of salt which prevents the bubbles from coalescing. These equations are taken from Sangani and Didwania 8 and Sangani et al. 6 As mentioned earlier, the velocity fluctuations induced by the mean gradient are much greater than the characteristic shear velocity based on the radius of the bubbles when Re is large. Thus, we will consider in Section III the simple case in which the bubbles have isotropic velocity fluctuations with no mean relative motion between the two phases, no imposed shear, and no viscosity. The properties of such bubble suspensions are functions of and the mean-squared velocity fluctuations characterized by the bubble-phase temperature T. Since this system is analogous to a molecular system in the absence of an imposed flow, we shall refer to this as the equi-librium state of a bubble suspension. We determine the equation of state, i.e., the relation between the bubble-phase pressure, , and T, using numerical simulations. The simulations are supplemented with a theory for dilute bubbly liquids. As shown recently by Yurkovetsky and Brady, 10 the equilibrium state properties of bubbly liquids can be determined from the configuration-dependent Hamiltonian of bubbly liquids using the standard statistical mechanical techniques. More specifically, we use the Hamiltonian to determine the pair probability distribution and hence the average properties of dilute bubbly liquids. The predictions of this theory are shown to be in excellent agreement with the simulation results. Since our primary interest is in developing a kinetic theory for sheared suspensions which are not in equilibrium and for which detailed pair probability distribution is not easy to determine, we also explore in this section how sensitive various properties are to the assumed expression for the pair probability density. It is shown that the simulation results are most consistent with an assumption that the dipoles induced by the bubbles are relatively independent of the distance between the bubbles.
In Section IV, we use Grad's moment method to develop a kinetic theory for sheared bubble suspensions. This theory accounts for the effect of the imposed shear on the velocity distribution of the bubbles in an approximate manner by considering only the second moments of the velocity distribution. Although this theory is similar to that for the gas-solid suspensions developed by Sangani et al., 4 the presence of lift force on the bubbles and the volume-fraction-dependence of the virtual mass lead to important modifications of the theory. The theoretical predictions are compared against the results of dynamic simulations of sheared suspensions in Section V and the agreement between the two is shown to be very good for a wide range of values of and Re. An important conclusion from this section is that the detailed hydrodynamic interactions are not critical in determining the rheology of sheared bubble suspensions. Rather, the role of potential flow interactions is only to set the average virtual mass of the bubbles and the viscous energy dissipation rate. These quantities may be determined as functions of from the simulations of the equilibrium state presented in Section III.
In Section VI, we investigate the influence of the finite Weber number of the bubbles on the results. The bubble deformation at a infinite Weber number tends to increase the added mass and drag coefficient of the bubbles. We assume an approximate expression for the ratio of the drag coefficient to the added mass that agrees with theoretical analyses of these quantities at moderate We and approaches the drag coefficient for a spherical cap bubble in the limit of very large We. The effects of bubble deformation on the dynamics of bubble-bubble collisions and the hydrodynamic interactions among the bubbles are neglected in this section. It will be seen that the effects of deformation on the drag result place an upper limit on the ratio of the bubble velocity variance to (␥a) 2 that can be achieved in practice. This limiting value is about 15. Nonetheless, the simulations with the nonlinear drag law confirm the possibility of achieving quite significant values of the bubble-phase pressure and viscosity.

II. A REVIEW OF BASIC GOVERNING EQUATIONS
In this section we review the basic equations governing the motion of bubbles in a liquid at Reӷ1 and Weϭ0. Further, we shall assume that the bubbles do not coalesce, but undergo elastic collisions. The coalescence of bubbles can be prevented in practice by the addition of surface-active impurities. Careful observations on the effect of surfaceactive impurities on the dynamics of pairs of bubbles have been made by  and Duineveld 14 who found that molar concentrations of surfactants as small as O (10 Ϫ4 ) are sufficient to prevent coalescence. Such low concentrations of the surfactants do not affect the potential flow approximation. At high surfactant concentrations, of course, the bubbles begin to behave like rigid spheres leading to breakdown of the potential flow approximation. Alternatively, the coalescence can also be prevented by the addition of salt as shown, for example, by Lessard and Ziemenski. 15 These investigators showed that a sufficient concentration of electrolyte in aqueous solution gives rise to short-range forces that prevent coalescence. Experimental observations of bubble collisions in salt solutions by Tsao and Koch 16 indicate that bubbles bounce with little loss of kinetic energy.
Let us consider N bubbles placed initially randomly within a unit cell of a periodic array with their centers at x ␣ and velocities w ␣ , ␣ϭ1, . . . ,N. The velocity of the bubble ␣ relative to the average suspension velocity The velocity of the liquid is similarly expressed as u(x,t)ϭ͗u͘ϩuЈ, where uЈ represents the disturbance flow induced by the bubbles moving with the relative velocity v ␣ . When the bubbles are spherical and the Reynolds number is large, the disturbance flow can be treated as irrotational. The velocity uЈ is then expressed as a gradient of a velocity potential satisfying with the boundary condition that n•ٌϭn•v ␣ on the surface S ␣ of the bubble ␣. Here, n is a unit outward normal vector on S ␣ . We shall use the method of multipole expansions to determine . As shown in Sangani and Didwania 8 and Sangani, Zhang, and Prosperetti, 17 the velocity potential can be determined accurately using a point-dipole approximation, where S 1 is a Green's function for Laplace equation in a periodic domain 18 and D ␣ is the induced dipole strength due to the presence of bubble ␣. The requirement that the average of uЈ over a unit cell of the periodic array be zero is satisfied by choosing 8 where is the volume fraction of the bubbles. Physically, G represents the back flow induced by the relative motion of the bubbles.
We shall treat the bubbles as massless and take the viscosity of the gas to be zero. Also, since we are interested in isolating the effect of the mean velocity gradient, we shall set the buoyancy force due to gravity to zero. The force balance on a representative bubble ␣ can be expressed in terms of its impulse defined by It is easy to show that the sum of the impulses of all the bubbles equals the momentum of the liquid induced by the relative motion of the bubbles. A bubble moving with a velocity that differs from the mean suspension velocity carries with it some liquid momentum, and therefore I ␣ may be thought of as the virtual momentum of bubble ␣.
The force balance on the bubble gives 8 where F p ␣ is the potential interaction force, F v ␣ is the viscous force, F ͗u͘ ␣ is the force due to mean flow, and F c ␣ is the force due to collision. The potential interaction force is evaluated from 19 where r represents the regular part of , i.e., the velocity potential minus the potential induced by the bubble itself. Therefore, the singular part of S 1 , i.e. 1/r, must be removed from S 1 before differentiating it three times for ␥ϭ␣. At large Reynolds numbers, the detailed short-time trajectories of the bubbles are determined primarily by the potential-flow interactions and collision forces. However, viscous forces play an essential role in controlling the kinetic energy of the suspension over long, O(Re Ϫ1 ), periods of time. To leading order, i.e., to O(Re Ϫ1 ), the viscous force can be evaluated by computing the rate of energy dissipation associated with the potential flow induced by the motion of the bubbles. Levich 20 used this observation to show that the drag on a single bubble is given by F v ϭϪ12av. His method was extended by Kok 11,12 to the case of many bubbles to show that the viscous force on individual bubbles can be determined from the relation F ␣ ϭϪ(1/2)ٌ v ␣E diss , E diss being the rate of viscous energy dissipation per unit cell. Alternatively, the viscous force on individual bubbles can also be determined by solving a Laplace equation for the viscous potential as shown in Sangani and Didwania. 8 As shown by these authors the two methods give identical results. Since the latter offers some computational advantages, we use it in the present study.
As in Sangani et al., 6 we shall use the following expression for evaluating the contribution of the mean flow to the force: where the derivatives of ͗u͘ are evaluated at xϭx ␣ , mϭ4a 3 /3 is the mass of the liquid having the same volume as the bubble, and D/Dtϭ‫ץ/ץ‬tϩ͗u͘•ٌ is the time derivative following the average motion of the suspension. While we do not have a rigorous proof justifying the use of ͑8͒, this relation is consistent with several known results for limiting and/or special cases as shown in the Appendix. These include ͑i͒ the case of small-amplitude oscillatory motion examined in Sangani, Zhang, and Prosperetti; 17 ͑ii͒ the linear extensional flow ͗u i ͘ϭe i j x j with e i j ϭe ji for which the flow is entirely irrotational; and ͑iii͒ the simple shear flow past a single spherical bubble with weak shear (␥aӶv ␣ ) for which the flow has a nonzero vorticity. 21 Finally, the expression is also consistent with that proposed by Auton, Hunt, and Prud'homme 22 for an isolated bubble if we The collision force is evaluated as described in Sangani and Didwania. 8 The collision is treated as an instantaneous event that preserves the kinetic energy and momentum of the liquid. This assumption is consistent with experimental observations and asymptotic analysis of bubbles bouncing in salt water at low We. 16 Upon collision with a bubble ␥, the impulse of bubble ␣ changes by F c k where k is the unit vector along x ␣ Ϫx ␥ and F c is the magnitude of the collision impulse determined from the relation 8 Here, w ␣ and w ␥ are the actual velocity of the bubbles ͓cf. ͑1͔͒ just before the collision and v ␣ and v ␥ are the velocities of the bubbles when each bubble is acted upon by a unit force along the line joining the center of the colliding bubbles and directed inwards toward the center of the bubble, the force on all the other bubbles in the suspension being zero.
Most of the results we shall present in this study were obtained by dynamic simulations in which the trajectories of the bubbles were evaluated using a modified Euler scheme described in Sangani and Didwania. 8 The simulations determined the microstructure parameters such as the radial distribution function, the velocity variance ͗v 2 ͘, and the dispersed-phase pressure tensor P i j ͑related to the bubblephase stress by i j ϭϪP i j ). As shown in Sangani and Didwania 19 and Bulthuis, Prosperetti, and Sangani, 23 the stress consists of three parts: kinetic, collisional, and potential interaction. These stress contributions are associated with the transfer of the ''virtual momentum'' associated with the bubbles across a surface by means of bubble translation, bubble-bubble collisions, and hydrodynamic interactions, respectively. They are defined by Here, nϭN/ is the number density of the bubbles, being the volume of the unit cell. The summation in ͑11͒ is over all the collisions occurring over time interval ⌬t.

III. PROPERTIES OF BUBBLY LIQUIDS IN THE ABSENCE OF SHEAR, GRAVITY, AND VISCOSITY
As mentioned in the Introduction, the root-mean-squared velocity fluctuations in the ignited state of a rapidly sheared bubble suspension are large compared with ␥a. Thus, it will be useful to consider the behavior of a bubble suspension in the absence of shear and viscous dissipation before developing the kinetic theory for the more practical situation of a sheared suspension of bubbles experiencing viscous forces. In the absence of viscous dissipation and shear work, the kinetic energy and momentum of the liquid are conserved and we will consider the case in which the liquid momentum and the mean velocity of the bubble-phase are both zero. In this case, the spatial structure and velocity distribution function of the system are isotropic and are determined entirely by the initial kinetic energy and the volume fraction of the bubbles. In analogy with the molecular systems, we shall refer to this as the equilibrium state of a bubble suspension. Our objective in this section will be to determine the equation of state, i.e., the relationship of the stress and kinetic energy to and the bubble-phase temperature Tϵ͗v 2 ͘/3. In addition, we shall be interested in comparing the microstructure of the bubble suspension with that of a hard-sphere molecular system, and in determining the leading order estimate of the rate of viscous energy dissipation when the liquid has small but finite viscosity. First, we will present results of numerical simulations for a range of volume fractions 0ϽϽ0.30. Then, in the second subsection, we will derive analytical results valid in the limit →0.

A. Numerical simulations
The spatial distribution of the bubbles is isotropic and can be characterized in terms of the radial distribution function g(r). Figure 1 gives results ͑diamonds͒ for the value of g at rӍ2a obtained by dynamic simulations of bubble suspensions. These values were computed by determining the number of pairs of bubbles with center-to-center separations in the range of 2a(1,1ϩ⌬) with the bin size 2⌬ϭ0.05. The results were obtained by averaging over several hundred a/T 1/2 time units in a simulation with Nϭ54. The time step for integrating the motion of bubbles was chosen such that on average a bubble underwent a collision in 50 to 100 time steps. Since many practical flows of bubble suspensions occur over low to moderate values of , we have computed the properties only for р0.3.
Also shown in Figure 1 are the results ͑pluses͒ for g(2a) obtained from a hard-sphere molecular dynamics code. We see that the values are very similar to those for the bubble suspensions. Thus, we conclude that the spatial distribution in the equilibrium state of bubble suspensions is very similar to that in a hard-sphere system. In particular, these results indicate that there is no significant clustering in this isotropic random motion of a bubble suspension. This may be contrasted with simulations of buoyancy-driven flows in bubble suspensions, which indicated that the potential flow interactions among the bubbles induced very significant clustering. 8,9 The solid line in Figure 1 corresponds to an estimate of ϭg(2a) obtained from the well-known Carnahan-Starling 24 approximation for hard-sphere molecular systems: The above equation is seen to give very good estimates of in both bubble suspensions and hard-sphere molecular systems. ͑It may be noted that the computed values are slightly lower than those predicted by the above expression because of the finite bin size used in the simulations.͒ We now present results for the bubble-phase stress. In the equilibrium state, the stress is an isotropic tensor characterized by a single scalar. The kinetic part of the stress is where (1/2)n͗I i v i ͘ϭ(3/2)n(m/2)C k ()T is the kinetic energy of the liquid per unit volume of the suspension. The virtual mass of an isolated bubble is m/2 and its impulse is I i Ӎmv i /2. Thus, C k will approach unity in a very dilute suspension in which bubble-bubble interactions are exceedingly rare. The results of dynamic simulations for various values of are shown in Figure 2. We see that the dependence of C k on is rather weak, and can be fitted adequately by a quadratic relation, Physically, (m/2)C k may be interpreted as a virtual mass of a bubble: it represents the energy required to increase (1/2) ϫ͗v 2 ͘ by a unit amount without applying any net force to the suspension. We shall refer to C k as the energy-added mass coefficient. It should be noted that this coefficient is very different from the momentum-added mass coefficient C a that previous investigators have examined. 17,25-27 C a relates the average impulse to the average velocity in the suspension and the coefficient of O() in its expressions reported in the literature ranges from 2.76 to 3.31 compared to approximately Ϫ0.35 found in the present case. When the net impulse and net velocity of the bubbles are changing with time, the test bubble experiences an effective medium whose density and relative acceleration differ from that of a pure liquid. This effective medium behavior alone contributes a correction of 3 to C a that is not present in the calculation of C k .
Since we have seen that the pair probability in a bubble suspension of volume fraction is very close to the value for the corresponding hard-sphere system, it is interesting to determine to what extent the stress in the bubble suspension can be related to the stress in an analogous hardsphere molecular gas. The kinetic stress in a gas of spherical molecules with mass m a can be written as 28 Comparing the above expression with ͑14͒, we see that the bubbles may be thought of as having a virtual mass m a equal to C k m/2. The kinetic stress in a bubble suspension is identical to that in a hard-sphere system with the same and kinetic energy. The collisional part of the pressure in a dense hard-sphere gas is given by 28 and therefore we define the collisional part in bubble suspensions in terms of a scalar C c given by with given by ͑13͒. The collisional stress would be identical to that in the hard-sphere system with the same kinetic energy and volume fraction if C c ϭC k . The values of C c determined from dynamic simulations are given in Figure 2.
Our simulations show that C c is in fact quite close to C k at the higher values of . For smaller , C c is somewhat smaller than C k with C c approaching approximately 0.97 as →0.
The relative importance of the potential and collisional parts of the stress is expressed in terms of a coefficient p defined by Both the potential and the collisional stresses are O( 2 ) for small and therefore p must approach an O(1) constant as →0. As shown in Figure 3 the potential stress contribution is rather small with p varying in the range 0.06-0.04 as varied from 0.05 to 0.3. In other words, the potential stress is only about 5% of the collisional stress and thus may be regarded as negligible. This very small potential stress in isotropically fluctuating bubble suspensions may be contrasted with the large potential stress calculated when there is a net relative motion between the bubbles and the liquid. In the latter case, the potential stress was larger than the kinetic and collisional stresses and it led to clustering of the bubbles and instability of the uniform state of the suspension. In bubble suspensions that contain a large velocity variance ͑due to shear or some other factor͒ as well as a relative motion of the two phases ͑such as that due to buoyancy͒, one may expect that the balance between the potential stress generated by the mean relative motion and the kinetic and collisional stress associated with the bubble velocity variance will determine the stability of the uniform state of the suspension. 10,[29][30][31] To develop the kinetic theory for sheared bubble suspensions, we shall need an estimate of the rate of viscous energy dissipation for a given variance in the velocity of bubbles. An isolated spherical bubble moving with velocity v through a nearly inviscid liquid experiences a viscous drag force equal to Ϫ12av. 32,20 We therefore express the rate of energy dissipation per unit volume of bubble suspension in terms of a scalar R diss defined by where R diss →1 as →0. The simulation results for R diss as a function of are given in Figure 2. It should be noted that although we evaluated the viscous force and dissipation at every time step to determine average dissipation rate, the trajectories of the bubbles were evaluated neglecting the viscous force on the bubbles so that the suspension microstructure depended only on the volume fraction of bubbles and the kinetic energy of the liquid. The results for R diss can be fitted by an approximate linear relation R diss is the viscous drag coefficient that relates the rate of viscous dissipation of energy in a suspension with no mean bubble velocity to the bubble velocity variance. It is different from the viscous drag coefficient C d defined, for example by Sangani and Didwania 8 and Sangani, Zhang, and Prosperetti 17 in terms of the mean viscous force acting on a suspension of bubbles with a mean motion relative to the suspension. The latter investigators found the coefficient of O() in C d to range from 2.11 to 2.28 compared with an approximate estimate of Ϫ0.16 in the present simulations for R diss . The effective medium due to mean relative motion, which contributes 2 to C d , makes no contribution to R diss . Not all properties of bubble suspensions can be interpreted in terms of a hard-sphere system with a volume fraction dependent virtual mass. For example, the collision frequency Ṅ c defined as the number of collisions per unit time per particle can be expressed in terms of T by means of For hard-sphere systems c f ϭ1. The results of numerical simulations for both bubble suspensions and hard-sphere suspensions are shown in Figure 4 where we see that while c f is indeed close to unity for the hard-sphere system the bubble suspensions yield c f equal to about 0.81 for the whole range of volume fractions. As the bubbles approach each other, their added mass increases and the relative velocity decreases leading thereby to lower collision rates. Also shown in Figure 4 is the coefficient for collision force determined in numerical simulations. This coefficient is related to the collision force by where k is a unit vector along x ␥ Ϫx ␣ . Note that m f c corresponds to the average of factor 2/(k•(v ␣ Ϫv ␥ )) in ͑9͒. According to the hard-sphere model, which we used earlier in describing the results for collision stress, f c would be simply equal to C k /2. The results of numerical simulations, however, show considerably greater magnitude of f c , and, in particular, as shown by the solid line in Figure 4, f c is approximately given approximately by 5/7C k . This additional factor of 10/7 arises from the added mass of the bubbles at collision being greater than that of isolated bubbles. Thus, apparently, the collision stress in bubble suspensions being approximately the same in the effective hard-sphere systems is a result of cancellation of two opposing effects of reduced collision frequency and collisional impact velocity on one hand and the increased added mass of bubbles at contact on the other.

B. Two-bubble calculations
In this section we determine the coefficients C k , C c , R diss , etc., in the limit →0. The calculations require the knowledge of the pair probability function P 2 (I 1 ,I 2 ,x 1 ,x 2 ) which represents the probability of finding two bubbles with their centers at x 1 and x 2 with impulses I 1 and I 2 , respectively. In the presence of dissipative effects, this pair distribution function must be determined from detailed trajectory calculations for pairs of interacting bubbles as has been done by Kumaran and Koch 33,34 and van Wijngaarden and Kapteyn. 35 These calculations, however, are generally tedious. Fortunately, this probability distribution can be easily determined based on a statistical mechanical argument as recently demonstrated by Yurkovetsky and Brady. 10 These investigators and Russo and Smereka 30 have shown that the dynamics of bubbles with potential flow interactions can be described in terms of a Hamiltonian H with the position and impulses of the bubbles treated as the generalized coordinates such that For bubbly liquids considered here H simply equals the kinetic energy of the liquid, i.e., Here, the configuration of bubbles C N consists of the positions and impulses of the N bubbles. The velocities of the bubbles are treated as dependent variables. For the special case of zero mean impulse, the probability density is given by the Boltzmann equation: 10,30 where A N and ␤ are given by The constant A N was determined by requiring that the N-particle probability must be equal to the product of N single-particle distributions when all the bubbles are well separated from each other and from the normalization condition ͐ P 1 (C 1 )dI 1 ϭn. The constant ␤ was determined from the calculation of ͗I 2 ͘ as explained later ͓cf. the discussion following ͑38͔͒. It should be noted that ␤ is related to the inverse of bubble-phase temperature based on impulse fluctuations and not the velocity fluctuations.
For two bubbles with center-to-center separation vector x 2 Ϫx 1 ϭaR, it is easy to show that the velocities of the bubbles are related to their impulses by where k i ϭR i /R is the unit vector along the line joining the centers of the bubbles. Here, and in all the subsequent calculations, we have neglected terms of O(R Ϫ9 ) and smaller. These higher order terms make a small contribution to the suspension properties even for touching bubbles, Rϭ2.
With the form of the pair probability distribution established, it is relatively straightforward to compute various equilibrium properties of dilute bubbly liquids. The radial distribution function is obtained by integrating the pair probability distribution in the impulse phase:

͑29͒
The integration in the above expression and in the subsequent calculations is facilitated by introducing In terms of these variables, the total kinetic energy of the liquid is given by where The integral on the right-hand-side of ͑29͒ is difficult to evaluate exactly for arbitrary values of R and therefore we shall evaluate it in the limit of large R by expanding the integrand in a Taylor series. Thus, for example, we write the part dependent on Ī as where we have neglected terms of O(R Ϫ9 ). Now substituting for the kinetic energy from ͑31͒ in ͑29͒, using dI 1 dI 2 ϭ8dĪdÎ, and making use of the Taylor series expansion ͑33͒ we obtain The integrals over the impulse subspace were evaluated with the help of the following formulas: Now substituting for A 2 from ͑27͒ into ͑34͒ and neglecting terms smaller than O(R Ϫ6 ) we obtain g͑aR͒ϭ1Ϫ 9 4R 6 . ͑37͒ According to this expression ϭg(2a)ϭ1Ϫ9/256ϭ0.965 at Rϭ2 indicating that the potential flow interactions cause only a slight depletion of the bubbles in the vicinity of a test bubble. This result also indicates that the added mass effect in the presence of random fluctuations in impulse introduces a small effective repulsive potential between pairs of bubbles as noted earlier by Yurkovetsky and Brady. 10 The numerical simulations described earlier ͑cf. Figure 1͒ showed g(2a) for bubble suspensions to be approximately the same as that for a hard-sphere system with the same volume fraction. The small 3.5% change in g(2a) predicted by the theory is smaller than the statistical errors in the simulations.
To determine C k , we need to estimate ͗I•v͘ and ͗v•v͘. Let be any dynamic variable ͑e.g. I•v) associated with the motion of bubbles. Then its average for dilute bubble suspensions can be estimated using

͑38͒
Here, 1 (I 1 ,x 1 ) is the value of for bubble 1 with impulse I 1 placed at x 1 and 1 (I 1 ,I 2 ,x 1 ,x 2 ) is the value for the same bubble given the presence of a second bubble with impulse I 2 at x 2 . Note that we have neglected the three-bubble interactions which will be unimportant in dilute suspensions. Also, the above expression for determining average quantities is valid provided that the difference 1 (I 1 ,I 2 ,x 1 ,x 2 )Ϫ 1 (I 1 ,x 1 ) integrated over the impulse subspace decays faster than ͉x 2 Ϫx 1 ͉ Ϫ3 . This is indeed the case in all the calculations presented in this section.
Since the impulses of the bubbles are independently specified, the second term on the right-hand-side of ͑38͒ is identically zero when 1 ϭ(I 1 ) 2 . Consequently, ͗I 2 ͘ is determined solely from the first term, and upon integration this yields the relation between ␤ and ͗I 2 ͘ given in ͑27͒.
The kinetic energy per unit volume of the suspension equals (n/2)͗I•v͘ and to evaluate it we substitute ϭI•v in ͑38͒:

͑39͒
The leading order term of O(R Ϫ3 ) in the above integral vanishes upon integration over the impulse space while the remainder simplifies to

͑41͒
Combining the above expressions with the definition of C k ͓cf. ͑14͔͒, we obtain This approximate estimate of the O() coefficient, Ϫ3/8ϭϪ0.375, is in very good agreement with the coefficient Ϫ0.35 obtained from numerical simulations ͓cf. ͑15͔͒.
The above results can also be used to estimate the viscous dissipation coefficient R diss . To determine viscous forces, we must solve for the viscous potential 8 which, with the dipole approximation, is given by The viscous dipoles are to be determined from the boundary condition n•ٌ v ϭϪ12D ␣ . Here, D ␣ is the induced dipole for the potential flow described earlier. The viscous force is given by F v ␣ ϭ4ad ␣ . On solving for d we obtain Since the velocity of a bubble is related to the dipoles by v ␣ ϭϪ(2D ␣ ϩD 3Ϫ␣ •ٌٌR Ϫ1 ), the rate of viscous energy dissipation per unit volume of the suspension is given by and upon combining with ͑40͒, ͑41͒, and ͑20͒ we obtain which is in very good agreement with the results of numerical simulations given by ͑21͒.
Next, we determine the collision frequency and the collision stress. The former is given by ͓cf. ͑22͔͒ where gϭv 1 Ϫv 2 is the relative velocity. It is easy to show that the relative velocity of two bubbles along the line joining their centers is related to their relative impulse by The above quantity and P(C 2 ) must be evaluated at Rϭ2 to determine the collision frequency ͑48͒. To keep the consistency in our calculations correct to O(R Ϫ6 ), however, we shall carry out integration first with arbitrary R and substitute Rϭ2 in the final result. Note also that the integration over g•kϾ0 is the same as the integration over Î•kϾ0. Now expanding the exponential in P 2 to O(R Ϫ6 ) as in the previous calculations, and using the following results: With Rϭ2, the above yields c f ϭ413/512ϭ0.807, in excellent agreement with the results of numerical simulations for the collision frequency shown in Figure 4. The first factor in ͑52͒ corresponds to ͗v•k͘/͗Î•k͘, which is the inverse of the factor by which the added mass of the bubbles changes due to interactions. The added mass of the bubbles at collision, i.e. at Rϭ2, is 32/23ϭ1.391 times that of isolated bubbles.
͑Note that this factor is slightly different from 10/7ϭ1.429 quoted earlier. The difference is due to terms of order smaller than R Ϫ6 retained in the derivation that gave the result 10/7.͒ This increased added mass, and consequently the decreased relative velocity, tends to lower the collision frequencies in bubble suspensions. The second factor in ͑52͒, which resulted from the integration in the mean impulse Ī-space, is the same as that obtained in the radial distribution function calculation. The effective repulsive potential induced by the potential flow interactions tends to deplete the pair probability density at contact and this too contributes to the decrease in the collision frequency. Finally, the third factor, which resulted from the integration in the Î-space, tends to increase the collision frequency. The terms of O(R Ϫ3 ) in P 2 (C 2 ) which did not contribute to the result for g(aR) make an important contribution to this last factor. Thus, a naive calculation, in which the pair probability is taken to be simply the product of two single-particle distributions with the radial distribution function, i.e., P 2 (I 1 ,I 2 ,x 1 ,x 2 ) ϭ g(aR)P 1 (I 1 ,x 1 )P 1 (I 2 ,x 2 ), would have given an incorrect estimate of 0.69 for c f . The collisional part of the pressure tensor is given by The collision force F c is equal to the relative velocity times the added mass at contact, i.e., Upon substituting for the collision force in ͑53͒, carrying out the integration, and comparing the result with ͑18͒ we find C c ϭ͓1Ϫ9/͑8R 6 ͔͒ Rϭ2 ϭ503/512ϭ0.982. ͑55͒ Once again, this result is in excellent agreement with the results of numerical simulations shown in Figure 2. Finally, the trace of the potential interaction part of the pressure tensor can be evaluated from

C. Sensitivity of the averaged properties to the pair probability distribution
We have seen that the average properties of dilute bubble suspensions obtained using the exact expression for the pair probability distribution for the equilibrium case are in very good agreement with the corresponding results obtained from the numerical simulations. Unfortunately, the equilibrium case considered in the present section is rather special and it is not easy to extend such exact analyses to nonequilibrium situations such as sheared suspensions. The theory in which the pair probability distribution can be factorized into one-particle distributions is much easier to develop for dense sheared suspensions and therefore it is of some interest to explore how various properties of bubbly liquids depend on the assumed form of pair probability. We have calculated the properties described above for dilute bubbly liquids for three simple forms of pair probability distributions: ͑i͒ P(C 2 )ϭn 2 f (v 1 ) f (v 2 ); ͑ii͒ P(C 2 )ϭn 2 f (I 1 ) ϫ f (I 2 ); and ͑iii͒ P(C 2 )ϭn 2 f (D 1 ) f (D 2 ). We shall refer to these as, respectively, the independent velocity, impulse, and dipole approximations. Thus, for example, we assume in ͑i͒ that the velocities of the two bubbles are independent of their separation vector aR and then determine their dipoles and impulses in terms of the velocities and R. The averaged properties are then determined by integrating the appropriate quantity in the velocity and R space. These approximations for various properties are summarized in Table I along with the results of the rigorous analysis presented in Section III B and the results of numerical simulations.
The evaluation of C k , R diss , and p does not require a specification of the form of the function f : the only restriction is the normalization condition, e.g., ͐ f (v)dvϭ1 for the case of an independent velocity approximation, or ͐ f (D)dDϭ1 for the independent dipole approximation. We see from Table I that all three properties are quite sensitive to the approximation made regarding the pair probability. For example, C k , which can be interpreted as an added mass that relates the energy of the suspension to the bubble velocity variance, is given by 1ϩ3/16 according to the independent velocity approximation, and 1Ϫ15/16 according to the independent impulse approximation. The corresponding results for the added mass coefficient for a collective acceleration of the bubbles C a are 1ϩ3.31 and 1ϩ2.76, respectively. 25,27,26,36,37 As mentioned earlier, the effective medium contributes 3 to C a but there is no corresponding contribution to C k . We shall see in the next section that the most important quantity controlling the rheology of bubble suspensions is the ratio C k /R diss . The independent dipole approximation yields C k /R diss ϭ1Ϫ/8 compared with 1Ϫ0.19 obtained by the numerical simulations. As seen in Table I, the independent dipole approximation gives the best estimate for this and the other properties of the bubble suspensions.
A clue as to why the dipole approximation gives superior estimates may be obtained by examining the ''head-on'' collision of two bubbles moving towards each other along the Ϯx 1 -axis. Consider two bubbles initially separated by a large distance moving towards each other with velocities Ϯv ϱ . The impulse and the dipole of the bubble moving in the positive x 1 -direction are given by m Ϫ1 I ϱ ϭϪD ϱ ϭv ϱ /2, and the total kinetic energy of the liquid is mv ϱ 2 ϭ4mD ϱ 2 . As the bubbles approach each other, their added mass or, equivalently, impulse increases in magnitude while the velocity decreases so that the total kinetic energy of the inviscid liquid remains constant. At the instant when the two bubbles come in contact, the velocity and impulse of the bubbles are related to their dipoles by v c ϭϪ͑7/4͒D c , I c ϭϪ͑5/4͒mD c . ͑58͒ Note that the added mass of the bubbles at contact is I c /v c ϭ5/7m compared with the added mass of m/2 for isolated bubbles -a result cited earlier in the discussion of numerical results for the collision frequency. Since the total kinetic energy remains unchanged for purely potential flow interactions, we have I c v c ϭI ϱ v ϱ ϭ2mD ϱ 2 from which we obtain Thus, the dipoles only change by about 4% during the head-on collision. The change in the dipoles for other relative orientations of the pair of bubbles is likewise expected to be small. Consequently, the pair probability distribution factorized in terms of independent one-bubble dipole distributions is more accurate in predicting the properties of bubbly liquids than factorizations based on velocity-or impulse-distributions.
In Table I we also gave results for the collision frequency and collision stress. To evaluate these quantities the detailed form of f is necessary. We assumed a Maxwellian form for f . Thus, for example, for the case of independent dipole approximation we chose,

͑60͒
where T Dc ϭ͗D c 2 ͘/3 is the temperature based on the dipole distribution at collision. Based on our calculations for the head-on collision of two bubbles we estimated T Dc to be given by which assumes that in the bulk ͗D 2 ͘ϭ͗v 2 ͘/4, a result based on an isolated bubble for which DϭϪv/2. Similarly, we chose T Ic ϭ(m 2 /4)(10/7)T for the independent impulse approximation and T c ϭ(7/10)T for the independent velocity approximation. With these forms for f we see from Table I that all three approximations gave reasonably good estimates of both the collision frequency and collision stress. In summary, the calculations presented in this section offer two ways to interpret the results of dynamic simulations. In the first, the bubbles are treated as hydrodynamically noninteracting particles with a virtual mass of C k (). This gives an estimate of the collision stress that is in reasonable agreement with the results of dynamic simulations. It, however, does not predict the collision frequency correctly. In the second, the bubbles are thought of as interacting hydrodynamically but with very little change in their dipoles. This gives reasonably good estimates of all the properties of the equilibrium state of the bubble suspension. In addition, we also described a third, more rigorous approach based on first determining the pair probability density and then determining various average properties. The last approach yielded quite accurate results for a dilute, equilibrium suspensions but its application to more general, dissipative bubbly flows and to higher bubble concentrations would require substantially more effort than the first two approaches.

IV. KINETIC THEORY FOR SHEARED SUSPENSIONS
In this section we develop a kinetic theory for a dense, sheared suspension of bubbles. We begin with the conservation equation for a dynamic variable associated with the motion of bubbles: where w j is the velocity of a representative bubble and is the rate of change of with time following the motion of the bubble. We shall assume that the mean velocity gradient is ␥ i j and write where x j is the center of the bubble and v i is the velocity relative to the mean flow. We shall restrict our attention to the case when there is no mean force acting on the bubbles and hence ͗v i ͘ϭ0. Our main interest will be to determine the velocity variance and the kinetic and collisional stresses.
Since the kinetic part of the pressure tensor is given by P i j k ϭn͗I i v j ͘, we apply the above conservation equation to i j ϭI i v j . Substituting for in ͑62͒ we obtain D Dt where D/Dt is the time derivative following the average motion of the suspension. Analogous to ͑6͒, the rate of change of i j following the motion of the bubbles can be decomposed into four parts due to potential-flow interac-tions, viscous forces, the mean flow, and bubble-bubble collisions. The results of the previous section suggest that the potential flow interaction force will have a minor effect on the behavior of a bubbly liquid with zero mean relative velocity and therefore we shall neglect the contribution from this force.
Next we will estimate the viscous contribution to the rate of change of the kinetic stress, ͗ i j ͘. It can be shown that ͗ i j ͘ and hence ͗ i j ͘ must be symmetric. The proof relies on the fact that the impulse of a bubble can be expressed as a product of a symmetric generalized added mass matrix and the velocity of all the bubbles in the suspension. 9 Thus, we shall estimate the contribution to i j from the forces exerted on the bubble due to viscous effects and the mean flow simply from their contributions to İ i by using ͗ i j ͘ ϭ ͗İ i v j ϩv i İ j ͘. In other words, we take ͗I i v j ͘ϭ͗v i İ j ͘. This approximation would be exact if the changes in the velocities and impulses of the bubbles had no effect on their spatial structure and therefore their added mass matrix. This approximation will become increasingly accurate with increase Re. In addition, we will use the large Re estimate of the viscous dissipation obtained in the previous section. Thus, the viscous contribution is expressed as The contribution due to the mean flow is written as The fluid acceleration term mD͗u i ͘/Dt in ͑6͒ does not contribute to the above term since it must be multiplied by ͗v i ͘ which vanishes in the present case of no mean relative motion. The above estimates of the contributions to i j from various forces acting on the bubbles do not require a detailed knowledge of the velocity distribution. A calculation of the contribution from the collisional force on the other hand requires a specification of the two-particle velocity distribution at contact, i.e., at Rϭ2.
In the theory of hard-sphere dense gases and granular materials it is usually assumed that the two-particle velocity distribution can be expressed as a product of single-particle distributions, i.e., P 2 ϭn 2 f (v 1 ) f (v 2 ). The collisional rate of change of any dynamical variable is then given by ͑see Jenkins and Richman 38 and Kremer and Rosa 39 for details͒ where v 1 and v 2 are the velocities of the colliding particles, gϭw 1 Ϫw 2 is the relative velocity, w ␣ ϭv ␣ ϩ␥•x ␣ is the actual velocity of the particle ␣, and and Ј are the values of for particles 1 just before and after the collision. The above expressions are appropriate to bubble suspensions provided that the velocities of the colliding bubbles are independent of each other. Our calculations in the previous section suggest that it is probably more appropriate to treat the dipoles of the bubbles as independent. Detailed calculations, however, show that the results obtained with the dipole approximation are essentially the same as those obtained from the above expressions provided that we treat the bubbles as rigid particles with independent velocities and mass (m/2)C k (). Thus, we shall evaluate the collisional contribution to ͗ i j ͘ using ͑67͒-͑69͒ by substituting Later we shall quote how much the results for S and Q would have been affected if we had used instead the independent dipole approximation for the colliding bubbles.
An exact solution for f would require a solution of the Boltzmann equation for the velocity distribution. However, this is cumbersome and therefore we shall use Grad's method 5 to obtain approximate solutions for the second moments of the velocity distribution. This method was shown to yield reasonably accurate results for gas-solid suspensions at finite Stokes numbers and finite inelasticity. 4 In this method, the velocity distribution is assumed to take a form where f M is the Maxwellian velocity distribution. It is easy to show that the constant a i j is related to the second moments of velocity by ͗v i v j ͘ϭT͑␦ i j ϩa i j ͒.

͑71͒
Since the bubble-phase temperature T equals one-third the velocity variance, we require that the trace a ii be zero.
Substituting ͑70͒ for f into the expressions for S and Q, taking ϭI i v j ϭ(m/2)C k v i v j , and evaluating the integrals yield ͑see Kremer and Rosa 39 ͒

͑73͒
where e i j ϭ(␥ i j ϩ␥ ji )/2 is the rate of strain tensor.
We mentioned earlier that the results based on the independent dipole approximation are essentially the same as ͑72͒ and ͑73͒. We now cite a few specific results. The coefficient of the T 3/2 term in the expression for S evaluated assuming that the dipoles of the colliding bubbles are independent of each other differs only by 7% from the above result based on the rigid particles with virtual mass of (m/2)C k . Similarly, Q i (v j ), which is same as the collision stress, was found to be essentially the same for the effective rigid particle and the independent dipole approximations.
We now restrict our attention to steady, homogeneous flows with ␥ kk ϭ0 for which ͑64͒ reduces to ͗Ṗ i j ͘ϭ0 since n͗P i j ͘ is independent of time and position. Substituting for the viscous, mean flow, and collisional contributions to Ṗ i j , and rearranging, we obtain where v ϭm/(24a) is the viscous relaxation time for an isolated bubble.
In the next section we shall compare the predictions of the kinetic theory with the results of dynamic simulations. The simulations with periodic boundary conditions are most conveniently carried out for the special case of simple shear for which ␥ i j ϭ␥␦ i1 ␦ j2 with x 1 -axis as the flow direction, we obtain the following four scalar equations for a i j from ͑74͒: where ϭ 24 5a 1/2 Solving ͑76͒-͑79͒ together with the condition that the trace a ii must be zero yields ϪSt Ϫ1 ͮ , ͑82͒ In addition, it is easy to show that a 13 ϭa 23 ϭ0. Finally, the condition that a ii ϭ0 yields the following cubic equation for : 3 Ϫͫ St Note that the above equation, which represents the fluctuation energy balance, determines the steady state bubblephase temperature. At this point, it is interesting to compare the theory for bubble suspensions presented here with the theory of gassolid suspensions at finite Stokes numbers due to Sangani et al. 4 In that the effective Stokes number St was defined as Stϭ␥m/(6aR diss ), m being the mass of the particles. The energy dissipation coefficient R diss for the gas-solid suspensions determined from the Stokes flow interactions is very different from that for bubble suspensions. In bubble suspensions there is an additional factor C k () which accounts for the variation of virtual mass of bubbles with . Another important difference is the contribution to i j due to mean flow ͓cf. ͑66͔͒. This can be decomposed into two terms: ͑i͒ due to spatial acceleration of the mean flow; and ͑ii͒ due to mean vorticity of the fluid. More specifically, the contribution to İ i from the mean flow is ͓cf. ͑6͔͒

͑86͒
where we have taken I i ϭm a v i , m a being the effective virtual mass of the bubbles. The last term on the extreme right side of the above expression, which is commonly referred to as the lift force term, was absent in the theory of gas-solid suspensions. Also, for gas-solid suspensions, the term inside the square brackets is replaced by the reaction force Ϫm͗u i ͘ϭϪmD͗u i ͘/DtϪmv j ␥ i j . When there is no mean relative motion, there is no contribution from the D͗u i ͘/Dt term, and consequently the term ␥ ki a k j ϩ␥ k j a ki in ͑74͒ replaces the ␥ ik a k j ϩ␥ jk a ki term in the gas-solid suspension theory. Comparing the cubic equation for T 1/2 derived here with that in Sangani et al., 4 we see that the first two terms in the two theories are identical. The other two terms have the same leading order behavior in the limit →0. Thus, the nondimensional temperatures in bubble suspensions and gas-solid suspensions at finite St are identical at low volume fractions. The stress components for the two suspensions, however, are unequal. In fact, it is easy to show that a 11 and a 22 for bubble suspensions correspond, respectively, to a 22 and a 11 in the gas-solid suspensions in the limit of →0. Thus the effect of lift force is to simply rotate the stress components by 90°in the plane of shear at small . When all three roots of the cubic equation ͑85͒ are real and positive, the theory predicts three steady states of which the state with the intermediate value of T can be shown to be unstable while the states with the lowest and highest T are stable and are referred to as the quenched and ignited states, respectively. ͑By stability here we mean the stability to small spatially homogeneous variations in the bubble-phase variance. Thus, for example, a small increase in the variance for the intermediate state will yield a greater increase in the energy input by shear than the increase in the viscous energy dissipation rate, and consequently, the system will move toward the ignited state. Likewise, a small decrease in the variance starting from the intermediate state will cause the system to become quenched. Whether these ignited and quenched states are also stable to small spatial variations in the bubble-phase volume fraction or temperature cannot be answered, of course, by this sort of reasoning͒. In the small limit, the analysis of Tsao and Koch 7 originally derived for gas-solid suspensions can be applied to bubble suspensions as well. These investigators showed that both states exist when StϾ24 1/2 and St 3 Ͻ1.5 and the final steady state depends on the magnitude of the initial velocity fluctuations in the suspension. On the other hand, the final state is the ignited state regardless of the initial conditions if StϾ24 1/2 and St 3 Ͼ1.5. Finally, only the quenched state exists for StϽ24 1/2 and St 3 Ͼ1.5. Since the relation ͑85͒ for T in bubble suspensions departs from that for gas-solid suspensions at finite , it is interesting to compare the behavior at finite . Let us denote by Re i the Reynolds number below which the ignited state does not exist and by Re q the Reynolds number above which the quenched state does not exist. The cubic equation ͑85͒ was solved numerically to determine these two critical Reynolds numbers as a function of and the results are shown in Figure 5. Multiple steady states are possible when Re i ϽReϽRe q . For ReϽRe i only the quenched state exists while for ReϾRe q only the ignited state exists. We see that these two curves meet at Ӎ0.014. Above this value of the cubic equation ͑85͒ permits only one real root and hence only one steady state.
For a given when Re is decreased from a high value for which the ignited state exists, the velocity variance of the bubbles will decrease smoothly until ReϭRe i below which there will be a sudden decrease in the variance because of the transition to the quenched state. Figure 6 shows the jump in the variance that will occur at ReϭRe i . Once again, at Ӎ0.014, the jump in the variance vanishes and hence T will be a smooth, continuous function of Re for greater than this critical value. This behavior of bubble suspensions is qualitatively similar to that in gas-solid suspensions where the critical was found to equal approximately 0.058. The present analysis shows that the multiple steady states described in Sangani et al. 6 occur only for very dilute suspensions.

V. COMPARISON WITH NUMERICAL SIMULATIONS
We now compare the predictions of the kinetic theory presented in the previous section against the results of numerical simulations for bubble suspensions subject to a simple shear flow ␥ i j ϭ␥␦ i1 ␦ j2 . Two kinds of suspensions were simulated: ͑i͒ bubble suspensions with full hydrodynamic interactions as described in Section II; and ͑ii͒ an effective hard-sphere model of a bubble suspension. In ͑ii͒ the bubbles were treated as if they have a virtual mass of (m/2)C k . The trajectories of the particles were evaluated using where F i c is the hard-sphere like collision force. Note that the detailed potential flow interactions are neglected in these simulations and C k and R diss are the same as that determined in Section II. It should be noted that the kinetic theory we developed in the previous section corresponds exactly to the situation described by the effective hard-sphere simulations. Thus, a detailed comparison of the results from ͑ii͒ with the theory will allow us to assess the validity of the Grad's approximation used in the kinetic theory while the comparison of ͑i͒ and ͑ii͒ will indicate the validity of replacing the detailed hydrodynamic interactions among the bubbles with volume-fraction-dependent virtual mass and drag coefficients.
All simulations were carried out with 54 bubbles ͑or particles͒, which were initially randomly placed inside a unit cell of a periodic array. The simulation with full hydrodynamic interactions was carried out typically over several thousand collisions. The number of collisions per bubble once the steady state was attained was thus typically in the range of 100-500.
In what follows we shall present results for a wide range of values of Re with the goal of assessing the kinetic theory. Obviously the assumptions made in evaluating the bubble trajectory, viz. small Weber number and large Re, will not be expected to apply for such a wide range of values of Re. We shall defer the question of determining the range of Re for which the model simulations carried out here are most likely to apply to Section VI. Figure 7 shows the velocity variance as a function of for Reϭ180. At such a high Reynolds number, the variance is much greater than ␥ 2 a 2 and the velocity distribution should resemble that of the equilibrium state studied in detail in Section II. Thus, we expect the Grad's approximation to be reasonably accurate. We see that the simulation results are generally in good agreement with the predictions of the theory with the maximum deviation of about 15% at ϭ0.05 and 0.3. We note that the velocity variance goes through a minimum around ϭ0.15. The viscous energy dissipation per unit volume of the suspension varies approxi- mately linearly with while the shear energy input, which roughly equals s ␥ 2 , approaches a constant value times T 1/2 as →0. At steady state the balance of the two requires that T/␥ 2 a 2 vary as (Re/) 2 at small . Here, s is the bubble-phase viscosity as explained in more detail below. For a fixed value of T, the contribution from the collision stress at given T increases s at rate that is faster than the increase in viscous dissipation R diss and, consequently, the velocity variance goes through a minimum.
The results for shear viscosity of the bubble-phase as a function of at Reϭ180 are shown in Figure 8. The shear viscosity is related to the bubble-phase stress by 12 ϭϪP 12 ϭ s ␥ 12 . ͑88͒ Since the calculation of the potential interaction stress is computationally intensive and since this stress was found to be roughly 5% of the collisional stress for the equilibrium case considered in Section III, we calculated the stress in the dynamic simulations by adding only the kinetic and collisional components of the stress. These stresses are evaluated in the theory from

͑90͒
As shown in Figure 8, the agreement between the theory and simulations is excellent when the results of dynamic simulations are plotted as s /T 1/2 . Thus, the 15% deviation in T between the theory and simulations noted earlier ͑cf. Figure 7͒ does not arise from errors in the expression for the shear viscosity in terms of T. At large Re, the bubble-phase viscosity is sensitive to the coefficient of the leading order terms in S and Q in ͑72͒ and ͑73͒, and a good agreement with the results for s indicates that the expressions for these terms based on a hard-sphere model are reasonably accurate.
The dotted line in Figure 8 corresponds to the shear viscosity based on the well-known kinetic theory of dense gases with the molecules of the gas being treated as having an effective mass (m/2)C k : s ϱ ϭ 8 5 1/2 aT 1/2 2 C kͫ 1ϩ We see that the calculated value of the shear viscosity is in reasonably good agreement with the above formula for Re→ϱ. The theory of dense gases considers small perturbations to the equilibrium state (Tӷ␥ 2 a 2 ) and therefore the above expression for s is relatively simple and only a function of T and . In contrast, the approximate theory we have developed here is more complete in the sense that it attempts to determine the viscosity even when ␥a is comparable to T 1/2 . Figure 9 shows the comparison between the theory and simulations for the bubble-phase pressure defined as onethird the trace P kk . We see once again an excellent agreement between the theory and simulations.
The results for the normal stress differences are shown in Figures 10-11. According to the theory, P 11 ϭ P 33 and P 22 Ͼ P 11 in the limit of small . The theory for finite predicts that P 22 Ϫ P 11 will change sign at Ӎ0.17. The lift force tends to increase the velocity fluctuations in the x 2 -direction and this leads to a greater kinetic stress for the 22 components. This results in the positive value of P 22 Ϫ P 11 at low . At higher volume fractions the significance of the collision stress increases and the collisions increase the 11 components of the stress because the imposed simple shear flow induces a relative motion of the bubbles along the x 1 -axis. The results of numerical simulations are seen to follow this general trend. It may be noted that this behavior of bubble suspension is different from that found in gas-solid suspensions where the lift force was absent and the mean flow effect was to induce greater fluctuations in the x 1 -direction. As a consequence, P 11 for the gas-solid suspensions was greater than P 22 and P 33 for the whole range of .
The calculations discussed above were carried out at large Re to test the accuracy of the expression ͑72͒ for S which could not be tested independently through the simulations of the kind described in Section III. At such large a Re, however, it is likely that the Weber number will also be large and the spherical bubble assumption cannot be justified. Thus, it is of practical importance to compare the theory and simulations for smaller values of Re. Figure 12 shows the velocity variance as a function of Re for ϭ0.3. The results based on hydrodynamically interacting bubbles are indicated by diamonds while those based on the hard-sphere model are indicated by pluses. We see that both are in very good agreement with each other and with the predictions of the theory. The same is true for the results for the shear viscosity shown in Figure 13. It is interesting to note that s /aT 1/2 remains approximately constant as Re is varied from 200 to 60. This constant value is in good agreement with ͑91͒ based on the kinetic theory of dense gases which predicts that the constant is equal to 0.326 for ϭ0.3. Similarly, the results for the bubble-phase pressure shown in Fig   In contrast to the above results, the agreement for the normal stress differences shown in Figures 15-16 is not as good. Our theory predicted the normal stress difference P 22 Ϫ P 11 to be negative while the simulations gave positive values. The stress difference, however, is less than 5% of the pressure and is therefore unimportant from a practical point of view. The results for the stress difference P 22 Ϫ P 33 on the other hand are larger in magnitude and seen to be in qualitative agreement with the predictions of the theory. Note also that the results calculated including the full hydrodynamic interactions are essentially the same as that obtained from the hard-sphere model for Re greater than about 60. The stress difference for smaller Re is greater for bubble suspensions than for the hard-sphere model.
The results for ϭ0.15 are shown in Figures 17-21. We see once again that the results for the velocity variance are in very good agreement with the theory for a wide range of Re. The bubble-phase shear viscosity scaled with aT 1/2 is seen to vary only by about 30% as Re is decreased from 200 to 40. The agreement between the theory and simulations for the normal stress difference is good for P 22 Ϫ P 33 but not for the other stress difference. The results for the hard-sphere model system are in excellent agreement with those for the bubble suspension. Thus, the observed discrepancy between the computed normal stress differences and the theory arises due to inaccuracies in the Grad's approximation.
Finally, the results for ϭ0.05 are shown in Figures  22-26. Although there appears to be a slight systematic difference in the velocity variance and shear viscosity for the FIG. 14. Bubble-phase pressure P scaled by T as a function of Re for ϭ0.3. The kinetic, collisional, and total pressures for the bubble suspensions are denoted, respectively, by squares, pluses, and diamonds. The corresponding results for the effective hard-sphere model are indicated by crosses, triangles, and stars. The lines represent the theoretical predictions. bubble suspension and the hard-sphere model system at the higher values of Re, the agreement of both with the theory is quite reasonable. Note also that, in contrast to the results at higher , the scaled viscosity s /aT 1/2 decreases substantially as Re is decreased from 200 to 40. The simple approximation based on kinetic theory of gases ͓cf. ͑91͔͒ does not exhibit this Reynolds number dependence and is not an adequate description for dilute bubble suspensions. On the other hand, the theory presented here which accounts for the change in the velocity distribution due to imposed shear adequately describes the rheology of bubble suspensions at small . Finally, we note that the agreement between the simulation results and the theory is very good for the normal stress differences. It may be noted that small theory predicts that P 11 ϭ P 33 while both our simulations and the finite theory show P 11 Ϫ P 33 to be significant. We also note that the results for ϭ0.05 show no evidence of an abrupt change in the behavior of bubble suspensions from the ignited state to quenched state found for very small . This is consistent with Figure 8 which shows that such a transition only occurs for Ͻ0.014.

VI. BUBBLES WITH NONLINEAR DRAG
In Sections II-V, we considered the ideal case in which the bubbles are spherical and they produce a fluid flow that may be described using the potential flow approximation. In addition, the drag coefficient for an isolated bubble was assumed to have a linear dependence on the bubble velocity, i.e., FϭϪ12va. While these approximations are reasonable for bubbles traveling close to their terminal velocity over a narrow range of bubble sizes, they would not hold for the wide range of bubble velocities predicted in the ignited state of a sheared bubble fluid. Thus, it is imperative that we assess the influence of nonlinear drag on the dynamics of sheared bubble suspensions to determine whether a state qualitatively similar to the ignited state predicted in the foregoing analysis can exist under physically realistic conditions. We have shown that hydrodynamic interactions play a modest role in the dynamics of sheared bubble suspensions at least in the regime for which the potential flow approximation is accurate. For this reason and because no theory is available to describe the hydrodynamic interactions in the more general case, we will neglect such interactions in this section and will consider only the nonlinear drag acting on the bubbles and bubble-bubble collisions.
Our knowledge of the drag on bubbles comes from extensive studies of the terminal velocity of bubbles in a variety of liquids by Haberman and Morton 40 and others. The Morton number, M ϭg 4 /s 3 , is a dimensionless number that depends only on the viscosity and density of the liquid, the interfacial tension, and the gravitational acceleration g.
Thus, this number can be used to categorize various liquids in terms of the type of buoyancy-driven bubble motion that will be obtained. The most important parameter leading to a large variation of the Morton number among various liquids is the viscosity. The Morton number can be as large as 10 5 in highly viscous oils and, at such large Morton numbers, the bubble is always deformed whenever inertia is important. However, water has a Morton number of order 10 Ϫ10 . In water and other low Morton number fluids, there exists a range of bubble sizes in which the Reynolds number, Re v ϭva/, is large but the Weber number We v ϭ v 2 a/s is O(1) or smaller.
A bubble rising in water is spherical at low to moderate Reynolds number and begins to take on an oblate spheroidal shape at We v Ϸ0.5 which corresponds to about 0.5 mm radius bubbles and an approximate Re v of 130. As the bubble deforms, its drag coefficient and virtual mass both increase because the deformed bubble carries along with it more of the liquid. The increase in the drag coefficient increases the rate at which the kinetic energy associated with the bubbles' motion is dissipated, while the increased added mass increases the total kinetic energy in the liquid for the same bubble velocity variance or bubble temperature. Since the drag coefficient increases more rapidly than the added mass, the net effect of the deformation effects at finite We will be to decrease the bubble temperature. Moore 41 and Lamb 42 derived theoretical predictions for the aspect ratio, drag coefficient, and added mass of spheroidal bubbles at low and moderate We v . These predictions have been shown to be in reasonable agreement with experimental measurements of the terminal velocity and aspect ratio for We v Ͻ1.7 by Duineveld. 43 At larger Weber numbers, the bubble shape loses its fore-aft symmetry and becomes unsteady and the drag on the bubble increases sharply. Eventually, at We v у30, the bubble assumes a steady spherical cap shape and the drag coefficient takes on a constant value, i.e., C D ϭ2F/(v 2 a 2 )Ϸ2. 6.
To obtain a rough estimate of the effects of the changes in drag and added mass induced by bubble deformation on the bubble velocity variance, we will perform simulations in which each bubble experiences a pseudo-steady drag force FϭϪ12Rva and has an added mass m a ϭ(m/2)C, so that the equation for the acceleration of bubble ␣ is where v ϭa 2 /(18) is the viscous relaxation time of a spherical bubble. This expression assumes that the coefficient for the lift force is the same as the added mass coefficient. Finally, although the bubbles are deformed at the higher Weber numbers, we will continue to treat them as elastic spheres when evaluating bubble-bubble collisions. We require a relationship for the ratio R/C over the full range of Weber numbers in order to model the behavior of a sheared bubble fluid with a wide distribution of bubble velocities.
Moore 41 provides expressions for R and C that are valid for We v Ͻ1.8. In the limit of small Weber number, these expressions yield R/Cϭ1ϩ0.038We v . Furthermore, C D /CϭRe v R/24C goes through a minimum at We v Ϸ1.3. However, Moore's expressions are not applicable at We v Ͼ1.8. We know from experimental observations of the terminal velocity that C D Ϸ2.6 in the limit We v →ϱ. There are no experimental or theoretical results for the added mass coefficient at larger We v , so we will simply assume that the added mass remains at the value calculated by Moore for We v ϭ1.8 for all higher Weber numbers. An empirical relationship for R/C that exhibits this large We v asymptote and approximates Moore's results for moderate We v is

͑93͒
In Figure 27, we present simulation results for the variance of the bubble velocity scaled with (␥a) 2 in suspensions with volume fractions of 0.15 ͑diamonds͒ and 0.3 ͑pluses͒. The ratio of the Weber and Reynolds numbers based on the shear rate is chosen such that We/Re 2 ϭ 2 /(a) Ϸ 1.4ϫ10 Ϫ5 , corresponding to bubbles with radii of 1 mm in water. It can be seen that the variance grows to a value 10 to 15 times larger than (␥a) 2 before the ratio ͗v 2 ͘/(␥a) 2 passes through a maximum at ReϷ70.
The kinetic theory result ͑74͒ for the second moments of the bubble velocity can be modified to include the effects of nonlinear drag, if the second term in ͑74͒ is replaced by the value of ͗v i v i ͘ derived using the new expressions ͑92͒ and ͑93͒ for the bubble acceleration. The integrals for the dissipation caused by the nonlinear drag must be performed numerically and so the equations for the moments cannot be written in an analytical form. However, they can readily be solved by Newton-Raphson iteration. The predictions of the kinetic theory with nonlinear drag are presented as dotted and solid lines in Figure 27. They approach the kinetic theory for spherical bubbles for ReϽ40 where nearly all the bubbles in the suspension have small values of We v . As the Reynolds number is increased, the theory exhibits a maximum value of ͗v 2 ͘/(␥a) 2 that is within about 20% of the maximum seen in the simulations. Both the theory and simulations indicate that the velocity variance is larger in the more concentrated bubble suspension as a result of the higher collision frequency.
Thus, we have seen that nonlinear drag has the effect of limiting the velocity variance that can be produced by shearing a bubble suspension to a value at which the Weber number based on the root-mean-square bubble velocity is of order one. In the parameter regime that can be easily obtained with millimeter-sized bubbles in low Morton number liquids, the nonlinearity of the drag associated with bubble deformation makes significant quantitative changes in the predictions of the theory. Nonetheless, the important qualitative prediction that shearing can induce a substantial variance of the bubble velocity and an associated bubble-phase stress remains valid in the presence of nonlinear drag which agrees with ͑8͒.