A RAISING AIR ALSO HAS A BUBBLE SYSTEM THAT WILL INFLATE
Rising gas bubbles are routinely used to agitate liquids of diverse composition ranging from pure water to molten metal, but bubble plumes themselves present a class of fluid dynamics problems that are complex enough to defy precise solution by any currently available means. Nevertheless, from an approximate standpoint, considerable progress has been made, first with analytical models , , ,  and  that employ pre-defined (Gaussian) distributions for bubble concentration and liquid velocity, and more recently with numerical simulations , , , , , , , , , , , , ,  and  that compute these quantities from prescribed governing equations. Single-phase simulations , , ,  and  (also called quasi single-phase simulations ) confine the rising bubbles to a predetermined axisymmetric plume, inside which an upward buoyant force is imposed on the surrounding liquid, with the latter regarded as a single continuous fluid. Two-phase simulations , , , , , , , ,  and  impose no such plume constraint and fall into two general categories. Eulerian–Lagrangian simulations , , ,  and  employ a Lagrangian viewpoint for the bubbles, tracking them as discrete particles, and an Eulerian viewpoint for the liquid, which is regarded as a continuous fluid. In contrast, Eulerian–Eulerian simulations , , ,  and  use a volume-averaged approach in which the liquid and gas are treated as distinct, but continuous, fluids. Balaji and Mazumdar  review single-phase models in particular, and Mazumdar and Guthrie  do likewise for single-phase and two-phase models in general. Jakobsen et al.  offer an extended review comparing the two different types of two-phase models.
Two-phase simulations may do an acceptable job of reproducing the near-field dynamics of bubble plumes, but they may also demand more resources than are feasible for computing the flow created by multiple plumes. In a practical setting, where computational fluid dynamics (CFD) might otherwise be used to investigate arrays of bubble diffusers for large-scale mixing applications, there arises the need for a plume model that can be easily incorporated into existing single-phase CFD codes. The present work accommodates this need with a simple approximation for bubble-induced vertical acceleration inside an idealized cylindrical plume, without solving additional equations of motion for the bubbles themselves.
Earlier single-phase models , , ,  and  have employed fixed cylindrical or conical plumes, within which all the bubbles remain as they rise toward the surface. In such models, the buoyant acceleration of the surrounding liquid is proportional only to the volume of undissolved gas inside the idealized plume, and Mazumdar and Guthrie  have demonstrated that approximations of this sort can render velocity predictions outside the plume that are as accurate as those of two-phase models. Although a conical plume is qualitatively more realistic, we propose a cylindrical model here because it is easier to implement for multiple plumes, and simplicity is a practical consideration when the model is to be used in designing arrays of two or more bubble diffusers to mix large volumes of liquid. For applications of the latter kind, in which velocity prediction is much more important outside the plume than inside, we demonstrate that a cylindrical column is an adequate representation for bubble plumes in unstratified liquids.
The present work implements a drift–flux model  for bubble plumes and a k–ϵ model  for turbulence in an existing single-phase CFD code  and . Empirical approximations are used for the plume radius  and the bubble slip velocity , but our numerical experiments indicate that these quantities have a relatively weak effect on the velocities outside the plume. Velocities computed with the proposed model are found to be most greatly affected by the total gas volume in the plume, and they are roughly proportional to the 2/5 power of the gas flow-rate from 200 to 22,000 cm3/s. Outside this range, however, we find that the exponent approaches 1/2 for much smaller flow-rates, and 1/3 for much larger flow-rates. The ability of the model to make useful far-field velocity predictions, with only the gas flow-rate, the liquid depth, and the container geometry as input, is demonstrated through comparison with experimental data for single and double bubble plumes in water.
- Cylindrical bubble plume
Bubbles released from an orifice quickly accelerate to their terminal (slip) velocities relative to the liquid , during which time they may also deform and break up into smaller bubbles . When this process is complete, we assume that the bubbles and the liquid are close enough to dynamic equilibrium that the drag force on the bubbles is balanced by the buoyant force on the liquid. The latter is given by Archimedes’ principle, and when the density of the gas in the bubbles is small compared with the density of the liquid in the bubble plume, the upward buoyant acceleration of the liquid in the plume is approximately
<img height=”33″ border=”0″ style=”vertical-align:bottom” width=”93″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si19.gif”>
where g is the constant acceleration due to gravity (9.81 m/s2) and αB is the local ratio of bubble volume to liquid volume. If the bubble volume ratio αB is also small, then Eq. (1) reduces further to
<img height=”12″ border=”0″ style=”vertical-align:bottom” width=”72″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si20.gif”>
Although αB can, in principle, take any value between zero and unity, we limit our consideration to situations in which αB≪1. In practice, values of αB greater than 0.1 progressively nullify the approximation given by Eq. (2), and they may also engender fluid dynamics too complex to be modeled adequately by a single-phase CFD code. In any case, if Eq. (1) were employed in favor of Eq. (2), it would produce larger buoyant accelerations (and correspondingly higher liquid velocities) than Eq. (2) with increasing αB.
Single-phase models , , ,  and  avoid computing the motions of individual bubbles by specifying the plume geometry in advance and by imposing simple approximations for αB. Dey Roy et al.  assume a cylindrical plume, and they assign a single value of 0.24 for the ratio of plume radius to container radius. Grevet et al.  and Mazumdar and Guthrie  assume conical plumes, with the cone dimensions inferred from laboratory experiments. Castillejos et al.  likewise employ conical plumes, but they also propose an empirical formula for the plume radius in terms of the gas flow-rate and the height above the bubble source. Inside the plume, single-phase model expressions for αB are given in terms of the gas flow-rate, either as empirical functions of position  or as drift–flux approximations ,  and  based on bubble residence times.
For simplicity, we idealize the plume as a right-circular cylindrical volume of liquid (Fig. 1) within which all the bubbles remain as they rise toward the liquid surface, and we derive an expression for αB based on the reasoning of Müller et al.  Let dVP be an infinitesimal cylindrical element of the plume, with radius rP and infinitesimal height dz, such that
<img height=”19″ border=”0″ style=”vertical-align:bottom” width=”100″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si21.gif”>
If dt is the time required for a bubble to traverse distance dz, then the ratio dz/dt is the local vertical velocity of the bubbles in the plume. This can be separated into two components, i.e.,
<img height=”35″ border=”0″ style=”vertical-align:bottom” width=”94″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si22.gif”>
with w the vertical velocity of the liquid, and w0 the slip velocity of the bubbles relative to the liquid. Defining Q as the volumetric flow-rate of gas bubbles through dVP, the bubble volume dVB enclosed within dVP at any time t is
<img height=”37″ border=”0″ style=”vertical-align:bottom” width=”166″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si23.gif”>
At a depth h below the liquid surface, the bubbles are compressed by the weight of the liquid, and the local gas-volume flow-rate through dVP is related to its value Q0 at atmospheric pressure by
<img height=”37″ border=”0″ style=”vertical-align:bottom” width=”88″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si24.gif”>
where h0 is the depth of liquid equivalent to one atmosphere of pressure (10.4 m for water). Combining , and with the definition for the local bubble-volume ratio in the plume, i.e.,
<img height=”37″ border=”0″ style=”vertical-align:bottom” width=”72″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si25.gif”>
the expression for the latter becomes
<img height=”39″ border=”0″ style=”vertical-align:bottom” width=”203″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si26.gif”>
Combining and , the buoyant acceleration of the liquid inside the plume is given approximately by
<img height=”39″ border=”0″ style=”vertical-align:bottom” width=”200″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si27.gif”>
Aside from the depth-correction factor, h0/(h+h0), Eq. (9) is the same as the drift–flux approximation introduced by Deb Roy et al.  Outside the plume, of course, the bubble-induced buoyant acceleration is identically zero.
<img class=”figure large” border=”0″ alt=”Full-size image (8 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr1.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr1.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr1.gif”>
Idealized bubble plume with radius rP in cylindrical container with radius R.
The quantities Q0, w, and h can be specified unequivocally in Eq. (9), because they are respectively determined by the bubble-generating device, the computed flow, and the height of the liquid surface. On the other hand, the plume radius rP and the bubble slip velocity w0 are free parameters in the proposed model, and they have to be determined empirically. If the flow velocity w were sufficiently large compared to the slip velocity w0, then the latter could be neglected in Eq. (9), but this does not hold for liquids accelerated from rest, for which w is initially zero. In any case, the assumption of a single fixed value for w0 has been widely used in single-phase and two-phase flow simulations , , , , , , , , , , , , ,  and , and we employ the same assumption here.
For bubbles ranging from 1 to 25 mm in diameter, the slip velocity varies from 15 to 35 cm/s . We choose the median (w0=25 cm/s) as our default value, and we demonstrate in Section 5 that the extreme values in this range (15 and 35 cm/s) produce little change in the liquid velocities computed with the proposed model. We show otherwise in Section 5 that the value chosen for rP does have a strong effect on the computed velocities, but only in the near field, at radial distances r≲2rP from the center of the plume. Since the plume radius must be determined empirically as a function of the gas flow-rate and the total depth, however, we postpone until Section 6 our discussion of specific formulas for rP.
- Liquid governing equations
The equations governing the motion of the liquid, both inside and outside the bubble plume, are taken to be the Reynolds-averaged continuity and momentum (Navier–Stokes) equations for incompressible flow, respectively given in vector form by
<img height=”14″ border=”0″ style=”vertical-align:bottom” width=”68″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si28.gif”>
<img height=”34″ border=”0″ style=”vertical-align:bottom” width=”215″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si29.gif”>
where ρ is the liquid density, uis the velocity, <img height=”11″ border=”0″ style=”vertical-align:bottom” width=”13″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si30.gif”> is the shear-stress tensor, gB is the buoyant acceleration created by the bubbles, and P is the pressure. The non-zero vertical component of gB is given inside the plume by Eq. (9), and the Reynolds-averaged tensor components τij of <img height=”11″ border=”0″ style=”vertical-align:bottom” width=”13″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si31.gif”> are represented as
<img height=”41″ border=”0″ style=”vertical-align:bottom” width=”171″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si32.gif”>
where the subscripts i,j may take values 1,2 or 3, corresponding to the three cartesian coordinate directions (x=x1, y=x2, and z=x3). Here the quantity ν indicates the kinematic eddy viscosity, rather than the molecular viscosity, and its value is computed from the k–ϵ turbulence model of Launder and Spaulding  via the standard formula,
<img height=”39″ border=”0″ style=”vertical-align:bottom” width=”74″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si33.gif”>
where k is the turbulence energy (per unit mass), ϵ is its dissipation rate, and the value of the empirical coefficient Cν is 0.09.
The turbulence model itself consists of two equations for the production, transport, and dissipation of k and ϵ; i.e.,
<img height=”38″ border=”0″ style=”vertical-align:bottom” width=”209″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si34.gif”>
<img height=”40″ border=”0″ style=”vertical-align:bottom” width=”273″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si35.gif”>
The product νΓ on the right sides of and is the production rate for turbulence energy, in which Γ is given by
<img height=”22″ border=”0″ style=”vertical-align:bottom” width=”432″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si36.gif”>
with the symbols u, v and w representing cartesian velocity components, and the subscriptsx,y and z indicating partial derivatives. The standard set of (non-dimensional) empirical coefficients for the k–ϵ turbulence model is C1=1.44, C2=1.92, σk=1.0, and σϵ=1.3.
Strictly speaking, the production rate for turbulence energy should include a direct contribution from the bubbles in addition to the conventional shear-production rate νΓ. In their Eulerian–Lagrangian model for bubble plumes, Johansen and Boysan  do include such a contribution, which arises from the relative motion between the bubbles and the surrounding liquid. For the proposed model, this contribution would be proportional to gBw0, making the total production rate νΓ+CBgBw0, where CB is a coefficient to be determined empirically. In practice, the additional production term itself presents no difficulty; but as far as we know, an adequate formula has not yet been developed for CB. Johansen and Boysan  present results for two gas flow-rates (1300 and 6100 cm3/s) with which they used two different values (0.1 and 0.7) for CB. They offer no general rule for specifying CB, but they do recommend further investigation of its possible relation to bubble size and shape , liquid surface tension, and turbulence length scales. An investigation of that sort is beyond the scope of our work here, but we do include in Section 7 a comparison of results obtained with CB=0 for a cylindrical plume, and CB=1 for a truncated conical plume. Otherwise, to avoid introducing an arbitrary free parameter in our model, we have excluded the direct production of turbulence by bubbles.
- Numerical flow solution
The proposed model for bubble plumes has been implemented in the MAC3D single-phase CFD code  and , which solves , , and in discrete (finite-volume) form. To achieve geometric flexibility in three dimensions, the equations are transformed from cartesian coordinates to general curvilinear coordinates  and  and discretized for structured curvilinear grids  created by a separate numerical grid generator . These grids may take familiar shapes like those obtained with conventional rectangular, cylindrical, and spherical coordinates; or they may take irregular shapes  that include rectangular, cylindrical, spherical, and distorted (boundary-fitted) components. In any case, the discrete flow variables are staggered in the grid cells, with normal (contravariant) velocity components defined on the cell faces, and scalar quantities (such as pressure) defined at the cell centers . The solution algorithm  is a semi-implicit time-marching scheme that uses second-order upwinding for advective terms , with the pressure gradient given by the iterative solution of a discrete Poisson equation  and  derived from and .
Concerning boundaries (Fig. 1), we idealize the upper (free) surface of the liquid as a frictionless lid (slip surface), where we impose null values for the vertical component of velocity and the vertical (z) derivatives of other variables (uz=νz=Pz=kz=ϵz≡0). The plume axis is treated as a slip line, along which we impose null values for the radial component of velocity and the radial (r) derivatives of other variables (wr=Pr=kr=ϵr≡0). Elsewhere, solid walls (container sides and bottom) are treated as frictional (no-slip) boundaries, along which standard wall functions  are used to define the values of k, ϵ, and shear stress.
In all the numerical simulations reported here, the flow was started from rest, with the plume acceleration gB activated at time t=0. Small initial values were assigned to k and ϵ, which made the initial eddy viscosity roughly equal to the molecular viscosity for water (10−6 m2/s). Except for the simulation of two interacting plumes in a large body of water ( Section 7), the code was allowed to run to steady state in each case, and the grids were sufficiently fine that further refinement (by a factor of two) changed the computed velocities by less than 3%.
- Numerical experiments
To investigate the behavior of the proposed model with respect to the plume radius rP, the slip velocity w0, and the gas flow-rate Q0, we used a hypothetical cylindrical tank with depth H=3.05 m, and radius R=15.25 m. The computational grid had a uniform spacing of 15.25 cm in the vertical and radial directions, out to radial position r=3.05 m, beyond which the radial spacing was gradually increased to 30.5 cm over an interval of 50 grid spaces. A hypothetical bubble source with a nominal gas flow-rate of 4720 cm3/s (at atmospheric pressure) was centered on the bottom of the tank (r=0, z=0).
Computed results are presented in non-dimensional form, using R and H to normalize the radial and vertical coordinates, r and z, respectively. Velocities are non-dimensionalized with the plume scaling velocity U0, which we define as
<img height=”46″ border=”0″ style=”vertical-align:bottom” width=”190″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si37.gif”>
Note that if U0 were indeed a characteristic velocity for bubble plumes, it would mean that the liquid velocity (predicted or observed) should scale with the 2/5 power of Q0 and the inverse square-root of H. The 2/5 exponent chosen here for Q0 lies roughly halfway between the extreme values of 0.33 and 0.48 observed experimentally  for bubble plumes in unbaffled containers filled with water. In any case, with the proposed 2/5 exponent for Q0, an inverse square-root of H is needed in Eq. (17) to give U0 in terms of Q0, H and g with proper units of length/time. Further discussion of the scaling relation between velocity and gas flow-rate is given in Section 8.
Fig. 2 and Fig. 3 show radial profiles (at mid-depth, z/H=0.5) of the non-dimensional vertical velocity w/U0, computed with three non-dimensional plume radii, rP/R=0.02, 0.03 and 0.04. Although the near-field profiles ( Fig. 2, r≲2rP) may look Gaussian, like those imposed in analytical models , , ,  and , their actual shapes lie somewhere between Gaussian and logarithmic, and they are the result of imposing the buoyant force given by Eq. (9) in the limited cylindrical region defined by r⩽rP. It is evident that the velocity near the idealized plume is quite sensitive to the value chosen for rP, with a ±33% variation in rP producing a ±30% spread in w. In the far field ( Fig. 3, r>2rP), the same variation in rP produces less than a ±10% change in w. In the numerical experiments discussed below, which demonstrate the separate influences of slip velocity and gas flow-rate, the plume radius was given the non-dimensional value, rP/R=0.03.
<img class=”figure large” border=”0″ alt=”Full-size image (7 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr2.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr2.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr2.gif”>
Near-field influence of plume radius rP on liquid velocities computed at mid-depth.
<img class=”figure large” border=”0″ alt=”Full-size image (8 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr3.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr3.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr3.gif”>
Far-field influence of plume radius rP on liquid velocities computed at mid-depth.
Fig. 4 and Fig. 5 show mid-depth radial profiles for w/U0, computed with three different values for the bubble slip velocity, w0=15, 25 and 35 cm/s. The latter span the range of values observed (in water) for individual rising bubbles ranging from 1 to 25 mm in diameter . In this case, a ±40% spread in w0 produces less than a ±5% variation in the computed velocity, both inside and outside the plume. This justifies the selection of the median observed value (w0=25 cm/s) as a default slip velocity for the proposed model.
<img class=”figure large” border=”0″ alt=”Full-size image (6 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr4.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr4.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr4.gif”>
Near-field influence of slip velocity w0 on liquid velocities computed at mid-depth.
<img class=”figure large” border=”0″ alt=”Full-size image (6 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr5.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr5.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr5.gif”>
Far-field influence of slip velocity w0 on liquid velocities computed at mid-depth.
Fig. 6 and Fig. 7 show mid-depth radial profiles for w/U0, computed with three different gas flow-rates: Q0=472, 4720 and 47,200 cm3/s. Recall that if Eq. (17) defines the characteristic velocity for bubble plumes, then w should scale precisely with the 2/5 power of Q0 and the inverse square-root of H. The results in Figs. 6 and 7 exhibit roughly a 2/5 power scaling with flow-rate over two orders of magnitude in Q0, both inside and outside the plume. As for the possible scaling of computed velocity with the inverse square-root of depth, we found (in other numerical experiments not presented here) that if both the gas flow-rate and the container geometry are held fixed, then w does indeed scale roughly as H−1/2. Fixed geometry means, however, that H/R and rP/H (or rP/R) must be the same from one simulation to the next. Otherwise, the scaling of velocity with depth may be quite different from inverse square-root.
<img class=”figure large” border=”0″ alt=”Full-size image (7 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr6.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr6.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr6.gif”>
Near-field influence of gas flow-rate Q0 on liquid velocities computed at mid-depth.
<img class=”figure large” border=”0″ alt=”Full-size image (7 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr7.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr7.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr7.gif”>
Far-field influence of gas flow-rate Q0 on liquid velocities computed at mid-depth.
- Bubble plume radius
Although the plume radius does not have to be known precisely for the proposed model to be useful outside the plume, it helps considerably if rP can be reliably estimated within, say, ±30%. This degree of uncertainty is small enough to limit the associated error in velocity to less than ±10% outside the plume (r>2rP). Thus, an empirical formula is needed, which roughly accounts for the variation of plume radius with the gas flow-rate and the vertical distance from the bubble source, so that a mean radius for the model plume can be specified in advance for any given flow-rate Q0 and depth H. For flow-rates Q0≲7000 cm3/s, Cederwall and Ditmars  propose the formula,
<img height=”31″ border=”0″ style=”vertical-align:bottom” width=”108″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si38.gif”>
with reference depth h0=10.4 m, and entrainment coefficient α taken as a function of Q0. Using values extracted by Cederwall and Ditmars  from data obtained by Kobus , we have developed an approximate expression for α, given by
<img height=”45″ border=”0″ style=”vertical-align:bottom” width=”142″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si39.gif”>
<img height=”46″ border=”0″ style=”vertical-align:bottom” width=”119″ alt=”” title=”” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-si40.gif”>
Note that the plume length-scale r0, first introduced by Castillejos et al. , was used in defining the scaling velocity U0 in Eq. (17). For application with the proposed model, however, Eq. (18) should be used with z=H/2 to obtain the mean radius of the bubble plume over total depth H.
- Comparison with experiment
Owing to the assumption of a stationary bubble column with a fixed radius, the proposed model can reproduce, at most, only the time-averaged flow created by real bubble plumes. This precludes simulation of phenomena such as wander , in which the plume gyrates about its (cylindrical) axis of symmetry. Even with this limitation, however, the model is still useful in computing velocity distributions needed for estimates of plume-induced mixing and transport . To demonstrate its capabilities in this regard, we compare model predictions with experimental data for single plumes at gas flow-rates of 205, 405 and 2500 cm3/s, all of which lie in the reported range of validity  for and . To demonstrate further applicability for multiple plumes, we also compare predicted and experimental results for the case of two interacting bubble plumes with gas flow-rates of 21,700 cm3/s each.
The first set of data comes from experiments done by Grevet et al.  in a water-filled cylindrical tank of depth H=0.6 m and radius R=0.3 m, into which air was injected through an orifice of diameter d=1.27 cm, centered on the bottom. The published velocity measurements consist of time-averaged velocity magnitude U=[u2+ν2+w2]1/2, obtained for two different flow-rates, Q0=205 and 405 cm3/s.
Fig. 8 compares our model predictions with the data for all six vertical positions (z/H=0.1, 0.3, 0.4, 0.56, 0.68 and 0.98) reported from the experiment . The computational grid had a uniform spacing of 1.5 cm in the radial and vertical directions. Using (to the nearest grid space) the values given by Eq. (18) with z=H/2=0.3 m, the model plume radii were specified as rP=7.5 cm for Q0=205 cm3/s, and rP=9.0 cm for Q0=405 cm3/s. Since two different flow-rates are involved, all velocities are presented in non-dimensional form U/U0, with U0 defined by Eq. (17).
<img class=”figure large” border=”0″ alt=”Full-size image (21 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr8.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr8.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr8.gif”>
Comparison of predicted and observed liquid velocities for cylindrical tank with depth H=0.6 m, radius R=0.3 m, and gas flow-rates Q0=205 and 405 cm3/s.
It is evident from Fig. 8 that the model is considerably less accurate in the vicinity of the idealized bubble column (r/R<0.3) than it is elsewhere. Otherwise, fair agreement is achieved between prediction and experiment, and the model generally reproduces the data outside the plume within a factor of two. The near agreement of the non-dimensional experimental results appears roughly to corroborate the proposed 2/5 power scaling of the velocity with Q0. With a factor-of-two ratio represented in the gas flow-rates, however, the difference between 1/2 and 1/3 power scaling is only about 12%. Any exponent in this range would probably render a fair correlation, given the scatter in the data.
The second set of data comes from experiments conducted by Castello-Branco and Schwerdtfeger  in a water-filled cylindrical tank of depth H=2.25 m and radius R=0.8 m, into which air was injected through an orifice of diameter d=5.0 mm, centered on the bottom. In this case, the published velocity measurements consist of radial profiles of time-averaged vertical velocity, taken at three vertical positions on opposite sides of the tank axis, for a single gas flow-rate, Q0=2500 cm3/s.
Fig. 9 compares our model predictions with the data for all three vertical positions reported from the experiments (z/H=0.089, 0.267 and 0.533). The computational grid had a uniform spacing of 2.5 cm in the radial and vertical directions. Using (to the nearest grid space) the value given by the Eq. (18) with z=H/2=0.6 m, the cylindrical plume radius was specified as rP=15 cm for Q0=2500 cm3/s. Since only one flow-rate is involved, all quantities are presented in dimensional form. Circular and triangular symbols indicate experimental data taken on the left and right sides of the tank axis,
respectively. Solid curves (—) indicate the proposed (simple) model with a cylindrical plume and standard k–ϵ turbulence model (CB=0). Dashed curves (- – – -) indicate the model (extended) with a truncated conical plume and direct bubble-generated turbulence (CB=1). In the latter case, the vertical variation of rP was approximated in “telescope” fashion, using values of 12.5 cm in the lower third, 15 cm in the middle third, and 17.5 cm in the upper third of the grid. This representation for rP exaggerates slightly the variation given by Eq. (18), but it approximates that variation to the nearest grid space, and it allows predicted results to be compared using the same computational grid for the cylindrical and conical plumes.
<img class=”figure large” border=”0″ alt=”Full-size image (16 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr9.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr9.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr9.gif”>
Comparison of predicted and observed liquid velocities for cylindrical tank with depth H=2.25 m, radius R=0.8 m, and gas flow-rates Q0=2500 cm3/s. The curves indicate model results obtained with (––) and without (—) a conical plume and direct bubble-generated turbulence. The symbols indicate experimental data taken to the left (•) and right (▴) of the container axis.
Outside the plume, the predicted results shown in Fig. 9 exhibit about the same degree of accuracy as those in Fig. 8. Unlike the previous case, however, the data of Castello-Branco and Schwerdtfeger  include a number of velocity measurements taken inside the plume itself. These reveal that the near-field distribution broadens, while the maximum velocity decreases, with vertical position (z/H). The computed near-field distribution also broadens, but the predicted maximum velocity increases with z/H for both the simple and extended versions of the model. The extended model reduces the predicted rate of increase somewhat, but not enough to achieve quantitative agreement inside the plume. Otherwise, there is little difference between the predictions of the two versions of the model outside the plume. Even if one allows for the data variation that arises here from plume wander, the model still fares poorly in the vicinity of the plume(r<0.2 m). The upshot is that the proposed model (simple or extended) is useful for predicting time-averaged velocities at radial distances greater than 2rP, but not for predicting velocities inside the bubble plume itself.
The final set of data comes from field tests conducted by Hornewer et al.  for coarse bubble diffusers in an abandoned, water-filled rock quarry, approximately 335 m long, 244 m wide, and 12 m deep. The experiments under consideration took place under destratified conditions with the diffusers placed near the middle of the quarry, 60 cm from the bottom. The case presented here involves two bubble diffusers, each with an air flow-rate of 21,700 cm3/s, separated by a horizontal distance of 4.6 m  and . Comparisons for other configurations with two diffusers and four diffusers can be found in a recent technical report .
For comparison with the velocity data obtained in the field test, the proposed model was executed (invoking symmetry) using the single-quadrant grid shown in the plan view in Fig. 10. Station 0 represents the left-side diffuser, station 1 the vertical axis of symmetry between the diffusers, and stations 2 and 3 locations respectively 2.3 and 13 m to the left of station 0. The computational grid extended well into the far field, to a radial distance of 50 m from the diffusers. The radial spacing was approximately 32 cm near the modeled diffuser, and the vertical spacing was uniformly 60 cm throughout the gird. The horizontal spacing was coarsened gradually with increasing distance from the diffuser.
<img class=”figure large” border=”0″ alt=”Full-size image (30 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr10.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr10.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr10.gif”>
Plan views of single-quadrant computational grid for two diffusers 4.57 m apart, showing left-side diffuser location (0) and experimental data stations (1, 2, 3).
The grid shown in Fig. 10 embodies the finest resolution that we tried for this particular problem, but it is not fine enough to remove grid errors entirely. We say this because factor-of-two refinement of cylindrical grids with comparable resolution may alter computed velocities by as much as 15%. In view of the model’s factor-of-two predictive capability, however, this level of possible grid error seems acceptable, given that a factor-of-two grid refinement demands at least a factor-of-eight increase in computer resources for three-dimensional flow simulations.Although the imposed air flow-rate (21,700 cm3/s) lies beyond the reported range of validity for Eq. (18), the latter nonetheless to specify a plume radius of 32 cm in the model input. Fig. 11 compares vertical profiles of velocity magnitude for stations 1, 2 and 3 after an elapsed time of one hour. Note that each of the experimental data represents a value that was time-averaged over a measurement interval of several minutes .
<img class=”figure large” border=”0″ alt=”Full-size image (15 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr11.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr11.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr11.gif”>
Vertical profiles of velocity after one hour for two diffusers 4.57 m apart.
The accuracy of prediction for this comparison is about the same, on average, as the previous comparisons with experiment. At stations 2 and 3 in particular, the predictions fall on either side of the data, roughly within a factor of two. The consistent underprediction of velocity along the vertical axis of symmetry at station 1, however, may be related to the model’s inability to account for plume wander with high flow-rate diffusers placed less than 5 m apart. In the absence of wander, this axis represents a line of minimum velocity with respect to neighbouring locations, and the computed magnitude increases sharply with horizontal distance from the axis. With wandering plumes, on the other hand, station 1 would no longer coincide with a fixed axis of symmetry, and the time-averaged velocities might be somewhat higher if plume gyration were included in the flow.
- Flow-rate scaling
In Section 5, Eq. (17), a reference velocity U0 was defined proportional to Q02/5, and it was pointed out that if U0 were indeed a scaling velocity for bubble plumes, then velocity in general should increase with the 2/5 power of the gas flow-rate. The flow-rate scaling exponent is of particular interest concerning the mixing times associated with plume-driven turbulent flow, because these may grow in inverse proportion to the liquid velocity .
For bubble plumes in water-filled laboratory models of metallurgical ladles, Mazumdar and Guthrie  report an exponent of 0.48 for flow-rates less than 500 cm3/s, and an exponent of 0.34 for higher flow-rates. In their empirical correlations for mixing times in general, however, Nakanishi et al.  propose an exponent of 0.4, while Asai et al.  propose an exponent of 0.33. Leitch and Banes , on the other hand, report a scaling exponent of 0.37 for velocities inside the plume itself. The exponents cited were all obtained for unbaffled containers, and their disparity suggests possible dependence on factors other than geometry. A survey of scaling relations for various configurations is given by Mazumdar and Guthrie .
With regard to the model proposed here, the scaling exponent is influenced by the relative magnitudes of the liquid velocity w and the bubble slip velocity w0 in Eq. (9). Note that the buoyant force that drives the flow in the plume is directly proportional to Q0 and inversely proportional to w+w0. The opposing shear force, which arises from non-uniformity in the local product of eddy viscosity and liquid strain rate, is dominated by a contribution proportional to w2 in the plume. If these forces are balanced at steady state, then the product w2(w+w0) becomes roughly proportional to Q0. Examining the limit as w0 approaches zero, one finds that w is proportional to Q01/3; but in the other limit as Q0 approaches zero, w is proportional to Q01/2. This asymptotic behaviour is qualitatively consistent with Mazumdar and Guthrie’s reported exponents of 0.48 for low flow-rates and 0.34 for high flow-rates. Whatever the case, given non-zero flow-rates and slip velocities, the scaling exponent must fall somewhere between 1/2 and 1/3.
To demonstrate the manner in which the predicted scaling exponent may vary with Q0, we used the proposed model to simulate plume-driven flow in a cylindrical container with 1 m depth and 2 m diameter. The computational grid had a uniform vertical spacing of 5 cm (20 grid spaces) and a radial spacing of 2.5 cm for the first 30 cm outward (12 grid spaces). Over the remaining 70 cm (18 grid spaces), the radial spacing was gradually increased from 2.5 to 7.9 cm. Flow simulations were executed for gas flow-rates from 1 to 1,000,000 cm3/s, with Q0 increased by a factor of 10 from one flow-rate to the next. The mean plume radius was given to the nearest grid space by Eq. (18), and the model was allowed to reach steady state for each case. Imposing the scaling relation UR∝Q0N, with UR denoting the maximum radial component of the return-flow velocity, we then extracted the value of the flow-rate scaling exponent N for each factor-of-10 increment in Q0. The results are plotted in Fig. 12, which indicates that N decreases gradually with increasing flow-rate, from a value approaching 1/2 for small Q0, to a value approaching 1/3 for large Q0. For the range of gas flow-rates (i.e., 200 to 22,000 cm3/s) considered in Section 7, the predicted exponent varies from 0.38 to 0.42, roughly approximating the proposed scaling exponent of 2/5 introduced in Eq. (17).
<img class=”figure large” border=”0″ alt=”Full-size image (4 K)” src=”http://origin-ars.els-cdn.com/content/image/1-s2.0-S0307904X9900030X-gr12.gif” data-thumbEID=”1-s2.0-S0307904X9900030X-gr12.sml” data-fullEID=”1-s2.0-S0307904X9900030X-gr12.gif”>
Predicted variation of flow-rate scaling exponent for single plume in cylindrical container with 1 m depth and 2 m diameter.
The drift–flux model that we have proposed for bubble plumes is similar to the cylindrical model employed by Deb Roy et al. , from which it differs mainly in its incorporation of a depth correction for the gas volume fraction. Owing to its simplicity, the model requires only a single additional term in the vertical (liquid) equation of motion, and it is easy to incorporate into existing single-phase CFD codes. Its utility is contingent upon neither the numerical scheme nor the turbulence model employed in the flow computation. The assumption of a fixed, empirically determined radius for the bubble column prevents the model from reproducing transient phenomena such a plume wander, and this reduces its accuracy very near the plume itself; but elsewhere the model generally reproduces observed (time-averaged) velocities within a factor of two, using as input only the gas flow-rate through the bubble source and the dimensions and geometry of the liquid container. This degree of accuracy may be inadequate for mixing calculations in the plume per se, but it is tentatively acceptable for the far-field mixing induced by single or multiple bubble columns in reservoirs.
In combination with a k–ϵ model for turbulence, the proposed model indicates that velocity increases in less-than-linear fashion with the gas flow-rate Q0. The predicted flow-rate exponent approaches 1/2 in the limit as Q0→0, and 1/3 in the limit as Q0→∞. This asymptotic behaviour is qualitatively consistent with some experiments , and our predicted exponent of 0.4±0.02, for flow-rates between 200–22,000 cm3/s, coincides with the median value in the range of exponents reported by others , ,  and . Since the predicted exponent depends on the relative magnitudes of the flow and slip velocities in the bubble plume, however, its value in general may vary somewhat with liquid depth and container geometry as well as flow-rate. Moreover, given the model’s factor-of-two limit in accuracy, its predictions for the flow-rate scaling exponent should be considered as rough estimates at best.
In its present form, the model is not applicable for stratified liquids in which there is a sharp density gradient between two weakly stratified layers. Our (unpublished) efforts to predict bubble-induced flow under such conditions have been unsuccessful , possibly because of the uniformly cylindrical model for the bubble column. Future work will investigate whether a more complex representation can be used successfully with a single-phase CFD code to compute the plume-driven flow under strongly stratified conditions.
To conclude all gases inside a bubble there can only be the center and the exit of the same amount of gas. There has to be the same volume of enter and exit.
Alvin L. Davis