## Abstract

We develop a model for the rheology of a three-phase suspension of bubbles and particles in a Newtonian liquid undergoing steady flow. We adopt an ‘effective-medium’ approach in which the bubbly liquid is treated as a continuous medium which suspends the particles. The resulting three-phase model combines separate two-phase models for bubble suspension rheology and particle suspension rheology, which are taken from the literature. The model is validated against new experimental data for three-phase suspensions of bubbles and spherical particles, collected in the low bubble capillary number regime. Good agreement is found across the experimental range of particle volume fraction (

## 1. Introduction

Multiphase suspensions of particles and/or bubbles in a continuous liquid phase are common in nature and industry; examples include magma, oil, concrete, foodstuffs, cosmetics, pharmaceuticals, biological fluids and nanofluids. Characterizing, modelling and controlling the flow of these suspensions requires a constitutive rheological model, encapsulating the viscosity of the suspension as a function of the properties of the suspending liquid, the volume fraction and properties of the suspended phase(s), and the flow conditions.

The rheology of two-phase suspensions (bubbles-in-liquid *or* particles-in-liquid, where the liquid is Newtonian) has been the subject of extensive experimental and theoretical research for more than a century. In recent years, significant advances have been made and two-phase constitutive equations are now available which have been validated against experimental data for a wide range of conditions (see [1] for a recent review). By contrast, considerably less research has been directed at understanding the rheology of three-phase suspensions (where bubbles *and* particles are suspended in a liquid) primarily owing to the complexity of the problem. Phan-Thien & Pham [2] present a theoretical treatment—discussed later in §3—which has been applied in studies of multiphase magma (e.g. [3,4]), but has not been experimentally validated. Experimental investigation of the rheology of three-phase suspensions appears to be confined to studies of bubble- and crystal-bearing magmas (e.g. [4,5]); these experiments and materials are complex and the resulting data are not well-suited to the validation of three-phase rheological models. Constraining three-phase rheology therefore remains an important, yet outstanding, problem in multiphase fluids research.

Here, we build on published two-phase constitutive equations to generate a three-phase model, by using an ‘effective-medium’ method in which the bubble suspension is treated as a continuous medium which suspends the particles; this carries the implicit assumption that the bubbles are small compared with the particles. We validate the model against new experimental data for three-phase suspensions of bubbles and spherical particles in the low-capillarity regime (in which flow is steady and bubble deformation is small).

## 2. Rheology of two-phase suspensions

The rheology of a strictly Newtonian fluid is completely described by its viscosity *μ*. The viscosity is the ratio of the deforming stress and associated strain-rate which, for rheometric flow, is given by *τ* is the shear stress and *η*_{r}, which is the apparent viscosity of the suspension at some strain-rate, normalized by the viscosity of the liquid phase

In the following subsections, we briefly review the constitutive equations for two-phase suspensions that provide the building blocks for the three-phase rheological model presented later in §3. The subscript ‘b’ refers to bubbles suspensions, the subscript ‘p’ refers to particle suspensions.

### (a) Bubble suspension rheology

When a bubble suspension flows, viscous stresses cause the bubbles to deform. If the flow is ‘steady’ the bubbles reach an equilibrium deformation, which is described by the bubble capillary number
*λ* is the bubble relaxation time [6–8]. The relaxation time describes the characteristic timescale over which the bubble adjusts towards a new equilibrium deformation in response to a change in the strain environment; it is given by
*a* is the bubble’s equivalent spherical radius and *Γ* is the liquid–gas surface tension. The flow is steady if the condition *t*≫*λ* [7].

For steady flow, the relative viscosity *η*_{rb} of a bubble suspension is given by Rust & Manga [8] and Mader *et al.* [1]
*η*_{b} is the apparent viscosity of the bubble suspension, and *η*_{r,0} and *Ca*. For non-dilute suspensions *η*_{r,0} and *et al.* [1]
*ϕ*_{b} is the bubble volume fraction. These expressions, which reduce in the dilute limit (as *η*_{r,0}=1+*ϕ*_{b}) and Mackenzie [11] *Ca* and decrease suspension viscosity at high *Ca*. Equation (2.4) is plotted in figure 1, which demonstrates that the transition between the asymptotic viscosity regions at low and high capillarity occurs over a fairly narrow range of *Ca*, centred on *Ca*≈1. Consequently, equations (2.5) and (2.6) can be used to calculate bubble suspension viscosity for all *Ca*, except over the narrow transitional region; we define an approximate upper bound to the low *Ca* region later in §4*c*.

Equations (2.2)–(2.6) are relevant for monodisperse bubble suspensions at low and moderate bubble volume fractions (*et al.* [1]. Bubble suspensions are visco-elastic even when dilute, and elastic behaviour becomes more pronounced as bubble volume fraction increases; visco-elastic rheology is neglected in this work because elastic behaviour is not manifest in steady flow [7].

### (b) Particle suspension rheology

Particle suspensions commonly show non-Newtonian behaviour when non-dilute, including shear-thinning (e.g. [12]), shear-thickening (e.g. [13,14]) and non-zero normal stress differences (e.g. [15, 16]). When shear-thinning behaviour is observed, the rheology of a particle suspension is often described using the model of Herschel & Bulkley [17]
*τ*_{0} is the yield stress, *K* is the consistency and *n* is the flow index (*n*<1 when the suspension is shear-thinning). The yield stress is non-zero only for highly concentrated suspensions, hence it is often neglected, reducing equation (2.7) to a power-law [18]
*η*_{p} is the apparent viscosity of the suspension. Although in common usage (e.g. [4,19,20]), this approach has the limitation that the consistency has fractional units of Pa s^{n} and is therefore not amenable to non-dimensionalization when *n*≠1; this issue is discussed in detail in Mader *et al.* [1] and Mueller *et al.* [20]. In this work, we address this limitation by introducing a characteristic timescale *t*_{c} of a shear-thinning suspension, against which the strain-rate can be non-dimensionalized, giving
*η*_{*} is a ‘reference viscosity’ of the suspension — i.e. the apparent viscosity at strain-rate *t*_{c} can be computed *a priori*. However, Mueller *et al.* [19,20] find empirically that the theoretical model of Maron & Pierce [21]
*ϕ*, *μ* and particle aspect ratio) when the consistency is identified with the viscosity; i.e. under the assumption *K*≡*η*. This is equivalent to finding that the characteristic timescale *t*_{c}=1 *s*, and making the identity *η*_{*}≡*η* in equations (2.9) and (2.10). This allows us to link these two equations while maintaining strict dimensional consistency, giving
*η*_{r,*} as the relative reference viscosity
*et al.* [19,20], Vona *et al.* [4] and Mader *et al.* [1], which included the pragmatic, but inexact, non-dimensionalization *K*_{r}=*K*/*μ*. Numerically, the values of *K*_{r} and *η*_{r,*} are identical when *t*_{c}=1 *s* (as is indicated empirically) so the results of those earlier studies can be transferred directly into this new framework.

The maximum packing fraction in equations (2.10) and (2.11) is a function of particle shape and roughness; Mader *et al.* [1] give the following equation for *ϕ*_{m}:
*r*_{p} is the particle aspect ratio. For smooth particles *ϕ*_{m1}=0.66 and *b*=1.08, and for rough particles *ϕ*_{m1}=0.55 and *b*=1.00; these values are empirically determined.

Mueller *et al.* [19] report that the flow index *n* for a particle suspension is a function of the particle volume fraction *ϕ*_{p} and the particle aspect ratio. They present a purely empirical relationship
*ϕ*_{p}/*ϕ*_{m}≤0.8.

## 3. A model for the rheology of three-phase suspensions

Equation (2.11) gives the relative viscosity *η*_{r,*} of a suspension of particles in a liquid with viscosity *μ*. If we suppose that the particles are instead suspended in a bubble suspension with viscosity *η*_{b} (i.e. we treat the bubble suspension as an ‘effective medium’) we obtain
*η*_{b}=*μ*(1−*ϕ*_{b})^{−1}; hence
*a*).

This effective medium method has been used elsewhere in rheological models. Of most relevance, Phan-Thien & Pham [2] use the approach to derive an equation for the viscosity of three-phase suspensions of bubbles and particles that is similar to the model we derive above, but contains a different expression for the particle suspension contribution: (1−*ϕ*_{p})^{−5/2}, which they derive using a differential method. Although they do consider a maximum packing fraction in some variants of their model, their implicit solutions reduce to exact equations only under restrictive conditions, e.g. *ϕ*_{p}∼*ϕ*_{m} or *ϕ*_{m}=1. Consequently, the treatment of the contribution of the particles to the suspension rheology in our formulation is a significant improvement over that of Phan-Thien & Pham [2].

### (a) Defining volume fractions in three-phase suspensions

Particle volume fraction and gas volume fraction are unambiguously defined for two-phase suspensions; however, care must be taken to define them appropriately for three-phase suspensions. In our model formulation above, we treat the bubble suspension as the effective medium, hence, the appropriate definitions are
*V* _{l}, *V* _{b} and *V* _{p} are the respective volumes of the liquid, bubble and particle phases.

For many three-phase applications, we are interested in characterizing how the rheology of a suspension of particles changes as bubbles are added to it (or, equivalently, as bubbles grow within it). For example, a magma that contains solid crystals may be bubble-free at depth, but become increasingly bubble-rich during ascent. In this case, it is more intuitive to define a particle volume fraction and bubble volume fraction as follows:
*ϕ**_{b} and *ϕ**_{p} that are known.

In the following sections, we work with both of these definitions, since equations (3.3) and (3.4) underpin the model formulation, while equations (3.5) and (3.6) are more natural for some applications of the model, and aid physical insight. A further volume fraction of interest is the fraction of the total volume that is made up of suspended bubbles and/or particles—the total suspended fraction:

## 4. Experiments

### (a) Samples

Three-phase samples were prepared by adding spherical glass beads (Potters Ballotini; density 2448 kg m^{−3}, size fraction 63–125 μm) to a sugar syrup (Tate & Lyle Golden Syrup; density 1438 kg m^{−3} and surface tension 0.08 N m^{−1} [7]) and aerating with a domestic electric whisk. The rheology of the pure syrup was determined individually for each sample batch and found to be strictly Newtonian; measured viscosities were in the range 55.68≤μ≤61.69 Pa s at 20°C (presented later in data table 1). Particle volume fraction was controlled by adding a known mass of beads to a known mass of syrup (typically equating to 100–150 ml) to prepare sample suites of similar initial (bubble-free) particle volume fraction *ϕ*_{p}=0.05, 0.1, 0.2, 0.3, 0.4 and 0.5. Bubble volume fraction was varied by adjusting the duration and speed of whisking, and suspension temperature. Errors in particle volume fraction and bubble volume fraction are ±3% and ±5%, respectively.

Of the resulting three-phase suspension, about 60 ml was used for rheometric analysis, and a small amount was imaged with a Zeiss SteREO V.8 stereomicroscope. With the remaining sample material, the bubble volume fraction *ϕ**_{b} was determined by measuring its weight and its volume in a 100 ml measuring cylinder. The bubble size distribution of each sample was determined using the image analysis software JMicroVision. A photomicrograph of a typical sample is presented in figure 2, along with its bubble size distribution.

### (b) Rheometry

Rheometric data were collected using a ThermoScientific Haake MARS II rheometer with Z40DIN concentric cylinder sensor geometry (rotor diameter 40.0 mm, cup diameter 43.4 mm, gap width 1.7 mm). A standard flow-curve determination consisted of a 20-step ‘up ramp’ of incrementally increasing shear stress *τ* up to a maximum value of 500 Pa (‘controlled-stress mode’), followed by a 20-step ‘down ramp’. At each stress step, the rheometer recorded the corresponding strain-rate

The densities of the particles and the suspending liquid are not well-matched in our experiments so settling must be considered; similarly, the bubbles are prone to buoyant rise (‘creaming’). The concentric cylinder sensor geometry was chosen because it is relatively insensitive to effects of settling and creaming compared with, say, a parallel plate geometry, because particles and bubbles move vertically past the sensor, rather than accumulating in a layer against it. From Stokes’ law, we can compute that the time required for an isolated particle or bubble in a dilute suspension to fall or rise the full length of the sensor is more than 120 h for the largest particle and around 1.5 h for the largest bubble; for a concentrated suspension, it is much longer because settling and creaming are hindered. Since total runtimes for our rheometric experiments are much shorter for each sample (the mean experiment duration was 5 min, maximum 12 min), the development of spatial gradients in the particle and bubble volume fractions is considered negligible.

### (c) Data analysis

Our rheometric experiments yield flow curves of applied shear stress *τ* against resultant shear strain-rate *Ca*′ for which the mismatch reaches 5%; i.e. for low capillarity

We determine the yield stress *τ*_{0}, consistency *K* and flow index *n* for each sample by fitting the Herschel–Bulkley model (equation (2.7)) to each filtered flow curve (figure 3) and determine errors using the bootstrapping method presented in the electronic supplementary material. For all samples, the yield stress is found to be either small and negative (which is unphysical) or positive, but within 2*σ* error of zero; hence, yield stress can be neglected and equation (2.7) can be expressed as the simple power law relationship given in equation (2.9). This is consistent with the experimental results of Mueller *et al.* [19], who found that yield stress is negligible for *η*_{*} and *n* in log-space using the relationship
*τ* and *b*, we assume *t*_{c}=1 *s*, based on previous experimental work [19,20]. The reference viscosity *η*_{*} that we determine is normalized by the viscosity of the syrup μ to give the relative reference viscosity *η*_{r,*} (equation (2.12)), hereafter referred to as the relative viscosity.

## 5. Results

Experimental data are presented in table 1. Results for relative viscosity are presented in figure 4, which plots *η*_{r,*}(*ϕ*_{p}), with datapoints coloured according to *ϕ*_{b}. Similarly, experimental results for the flow index are presented in figure 5, which plots *n*(*ϕ*_{p}), with datapoints coloured according to *ϕ*_{b}.

It is useful at this stage to confirm that the data for the two-phase end members (bubble-free particle suspension and particle-free bubble suspension) adhere to the relevant two-phase constitutive equations. The solid curve in figure 4 represents the best fit of equation (2.11) (in which the only free parameter is *ϕ*_{m}) to data for bubble-free particle suspensions (*ϕ*_{b}=0). We find an excellent fit for *ϕ*_{m}=0.593, with *R*^{2}=1.00; this value of the maximum packing fraction is slightly lower than the value of *ϕ*_{m}=0.633 quoted by Mueller *et al.* [19,20] for suspensions of monodisperse spheres, but we note that the 2*σ* errors of the two estimates overlap. Both of these values are lower than the value of *ϕ*_{m}=0.66 calculated from equation (2.13) with *r*_{p}=1. We propose that equation (2.13) overestimates *ϕ*_{m} for nearly spherical particles because it is based on a fit to data that assumes a Gaussian relationship between *ϕ*_{m} and *r*_{p}. Studies of non-sheared particle packs have reported that the maximum packing fraction is actually highest for slightly non-spherical aspect ratios [22,23]; consequently, the *ϕ*_{m}(*r*_{p}) curve dips around *r*_{p}=1 leading to the overestimate given by equation (2.13).

Figure 6 plots two-phase data for the relative viscosity of particle-free bubble suspensions, i.e. *η*_{r,*}(*ϕ*_{b}) with *ϕ*_{p}=0. The solid curve is equation (2.5), which gives a good fit to data, with *R*^{2}=0.87. Figures 4 and 6, therefore, demonstrate the validity of the two-phase constitutive models (equations (2.5) and (2.11)) on which we build our three-phase model.

## 6. Discussion

### (a) Reference viscosity

In figure 4, all suspensions that contain bubbles have a higher reference viscosity than a two-phase particle suspension with the same particle volume fraction. Conceptually, this is equivalent to saying that, for a given particle suspension, the viscosity increases if some of the suspending liquid is replaced with bubbles. This is intuitive, because bubbles in the low-capillarity regime increase suspension viscosity.

Recasting our data in terms of *ϕ**_{b} and *ϕ**_{p} yields figure 7. It is evident from this plot that the effect of *adding* bubbles to a particle suspension (or growing bubbles in a particle suspension) depends upon the initial particle volume fraction *ϕ**_{p}. For dilute particle suspensions (

These relationships are more clearly demonstrated by figure 8. The plot shows that the rheology of a particle suspension becomes increasingly sensitive to the addition of a small volume fraction of bubbles as its particle volume fraction approaches the maximum packing fraction. The physical explanation for this behaviour is straightforward and relies on two competing processes. As discussed above, the addition of low-capillarity (i.e. spherical) bubbles to a fluid increases its viscosity. For dilute particle suspensions, this is the dominant trend, hence our data show an increase in reference viscosity with increasing bubble content for *ϕ*_{p}, reduces the impact that particle–particle interactions have on suspension rheology and reduces suspension viscosity. This process dominates for concentrated particle suspensions; indeed, because the Maron–Pierce relationship is a power law, the higher the initial particle volume fraction, the greater the impact the same dilution with bubbles will have.

It is also clear from these figures that the data agree well with the three-phase model that we propose (equation (3.2)). The model predicts the relative viscosity to within ±20% for all but two of the samples, and within ±10% for the majority (figure 9). Furthermore, there is no strong systematic trend relating the discrepancy to either *ϕ**_{b} or *ϕ**_{p}. We note that the close agreement between model and data is found despite the fact that some samples violate the model assumption that bubbles are small compared with particles (§3). In all samples, the majority (by number) of the bubbles are smaller than the particles and so it makes sense to choose the bubble suspension as the effective medium. However, the average bubble radius is between 0.20 and 0.97 times the average particle radius and the bubble volume-mean-radius (§4*c*) is between 0.57 and 2.4 times the particle volume-mean-radius indicating that the ‘typical’ bubble in many samples is comparable in size to the particles.

The largest discrepancy between our model and data occurs when both *ϕ**_{p} and *ϕ**_{b} are large: our model tends to underpredict the reference viscosity in such samples. This discrepancy is probably caused by bubble–particle interactions, which are not captured in our simple approach of combining two-phase equations. Further work is needed to formulate a model that captures such interactions.

For a given bubble and particle volume fraction, it is useful to know whether the addition of further bubbles (or the growth of existing bubbles) will result in an increase or decrease in viscosity. An analysis of this scenario is presented in the electronic supplementary material.

### (b) Flow index

Reference viscosity provides only a partial description of the rheology of a shear-thinning suspension; for a practical rheological model, the flow index *n* is also required (equation (2.9)). Figure 10 re-presents the flow index data shown in figure 5, plotting them against *ϕ**_{b} and *ϕ**_{p}/*ϕ*_{m}.

Figure 10*a* shows clearly that shear thinning is observed for all suspensions, even those containing only bubbles. Our data indicate a linear relationship between *n* and *ϕ*_{b} for particle-free suspensions, in the low-capillarity regime, such that shear thinning becomes more pronounced as bubble volume fraction increases:
*ϕ*_{p}=0.) Bubble suspensions are known to be strongly shear thinning in the range *a*), which makes this result surprising. One potential explanation is that, despite our data filtering methodology (§4*c*), there are a small number of very large bubbles in the samples, larger than those measured in our sample images, which have a correspondingly long relaxation time and, as a consequence, have a capillary number in the transitional regime. An alternative explanation is that the current model for bubble suspension rheology (equation (2.4)) is inadequate for non-dilute suspensions. That model can be derived from the theoretical treatment of Frankel & Acrivos [24] and Llewellin *et al.* [25], which is analytically exact in the limit of a dilute suspension (in which bubble–bubble interactions can be neglected) and in the limit of small bubble deformations. It is possible that bubble–bubble interactions in non-dilute samples act to introduce shear thinning—this would be consistent with our finding that shear thinning becomes more pronounced as bubble volume fraction increases. Further experimental work would be required to underpin a more detailed investigation of this phenomenon.

Superimposed on the decrease in flow index due to increasing bubble volume fraction is the effect of increasing particle volume fraction, shown most clearly in figure 10*b*. This relationship is nonlinear and appears to follow the empirical model proposed by Mueller *et al.* [19] (equation (2.14) with *r*_{p}=1 for spherical particles). For bubble-free suspensions—comparable to those investigated by Mueller—there is excellent agreement between data and model for *ϕ*_{p}/*ϕ*_{m}<0.8, which is consistent with the limits of applicability given by Mueller *et al.* [19].

The combined effect of bubbles and particles on the flow index appears to be a simple superposition of these two effects: the flow index of a pure fluid is 1 and is reduced by some amount dependent on bubble volume fraction, and again by some amount dependent on particle volume fraction. Consequently, we propose the following purely empirical model for flow index for suspensions of spherical particles:
*b* and indicate that the model is valid for all samples with a total suspended fraction

For suspensions of non-spherical particles, the aspect ratio may be included in equation (6.2) in a manner analogous to equation (2.14). However, further experiments would be required to validate this extension.

## 7. Conclusion

Our results demonstrate that the proposed three-phase rheological model (equation (3.2)) based on an effective medium approach is in close agreement with experimental data over the full range

Mueller *et al.* [20] demonstrate that the Maron–Pierce relationship (equation (2.11)) is valid for suspensions of non-spherical particles when *ϕ*_{m} is calculated as a function of particle shape (equation (2.13)). Adopting this methodology for equation (3.2) broadens its applicability to natural systems, in which particles are rarely spherical.

The low-capillarity assumption can also be relaxed. Substituting equation (2.6) for equation (2.5) in the formulation of equation (3.2) yields a three-phase model suitable for high-capillarity flows. This version of the model predicts that the addition of bubbles will always reduce the reference viscosity of a three-phase suspension, even when particle concentration is low. At high particle concentrations, the reduction in viscosity is much more dramatic than in the low capillarity case.

Relaxing the assumption of steady flow is more challenging because, while the rheology of bubble suspensions in unsteady flow is well known [7,9], there is no adequate model for particle suspensions in unsteady flow.

The model developed in this study assumes that bubbles are small compared with particles, although the experimental data demonstrate that the model remains valid when bubbles and particles are comparable in size. Further experiments are required, however, to determine the rheology of suspensions in which bubbles are large compared with particles.

The proposed model for the flow index of three-phase suspensions (equation (6.2)) also provides good agreement with experimental data for which

## Data accessibility

The full rheological data can be found online at the Earth Science Academic Archive, doi:10.5285/e729d509-a616-4787-bb41-44c1169685b6.

## Funding statement

J.M.T. acknowledges support from NERC Studentship NE/K500999/1; S.P.M. acknowledges support from NERC Research Fellowship NE/G014426/1; we would also like to acknowledge support from RCUK in making this article open access.

## Acknowledgements

We would like to thank two anonymous reviewers for their helpful comments.

- Received July 22, 2014.
- Accepted October 28, 2014.

© 2014 The Authors. Published by the Royal Society under the terms of the Creative Commons Attribution License http://creativecommons.org/licenses/by/4.0/, which permits unrestricted use, provided the original author and source are credited.