Basic Concepts

Bubbles are Lagrangian spheres that can access continuous positions across the fluid lattice. They are used to model two-phase flow, gas transport processes, gas hold-up, and gasified power draws.

Bubbles enter the system along the children geometry of the bubble parent with a user-defined injection velocity. Bubble positions and velocities evolve according to Newton’s Second Law using a Verlet algorithm advanced at the simulation timestep. As discussed in Section Bubble Forces below, the forces included in the velocity update are defined by the user.

Upon entering the simulation domain, each bubble is assigned an ID, a birth timestamp, a diameter (see Section Initial Bubble Size Distribution below) and an origin ID. The birth timestamp identifies the time-step at which each bubble enters the system. This value is particularly useful in predicting bubble residence time distributions and bubble mean-age. The origin ID describes where a bubble entered the system . This value is useful in predicting how bubbles from various sources blend and transfer mass to/from the fluid. It is also useful in predicting how bubbles with different properties (e.g. density and diameter) are affected by fluid motion (assuming different bubble properties are assigned to the different bubble origins).

Bubble break-up and coalescence is discussed in Section Bubble Breakup and Coalescence. Breakup typically occurs on the level of individual bubbles and is informed by the physical properties of the bubble and the kinematic properties of the fluid. Coalescence typically occurs at the level of bubble pairs, and is informed by the kinematic properties of the two bubbles.

Mass transfer between bubbles and the surrounding fluid are discussed in Section Bubble Mass Transfer Coefficient and Bubble/Scalar Coupling. The mass transfer coefficient of each bubble is a user-defined function governed by local fluid properties. This bubble-specific mass transfer rate multiplied by the local bubble/fluid concentration difference defines the mass transfer rate between bubbles and the surrounding fluid. Mass leaving/entering the bubble causes the bubble to shrink/grow in a way that conserves total species mass.


Bubble output files are printed at the same frequency as output slices.



Name of the bubble set


Turns the bubble parent on/off in the solver input file

Injection Geometry

Two options are available

Inject at surface

Bubbles will be injected along the surface of the children geometry.

Inject in volume

Bubbles will be injected within the volume enclosed by the children geometry.

Injection Option

Six different injection options are available


Remove all bubbles that enter the geometry


Maintain a minimum user-defined bubble population inside the children geometry over a period


Dump a user-defined number of bubbles into the children geometry at a given time


Continuously feed bubbles into the children geometry at a user-defined number flow rate over a period

Volume Feed

Continuously feed bubbles into the children geometry at a user-defined volume flow rate over a period

Injection Number Scale

Scale factor applied to the number of bubbles entering the system, [-]

For large-scale systems with high volume feed rates, which may have billions of bubbles, the computational requirements necessary to monitor every bubble may exceed available memory. To make such simulations tractable, users may specify an Injection Number Scale that reduces the number of explicitly tracked bubbles by this factor. The simulated volume of the bubbles, along with hold-up, mass transfer, and two-way coupling effects, are then automatically multiplied by this factor inside the code. To maximize simulation fidelity, this number scale should be set to the smallest number necessary.

Resolved Feed Rate

Actual number of bubbles that will be fed into the system, [#/s]

The resolved feed rate describes the rate at which bubbles enter the domain through the children geometry. When the Injection Number Scale is set to 1, the resolved feed rate is equal to the user-defined volume feed rate divided by the mean bubble volume. Otherwise, the resolved feed rate will be smaller by a factor equal to the Injection Number Scale. Note that bubbles are injected uniformly across the entire child(ren) geometry surface area (or volume, if injecting in the child volume).

Start Time

Time to begin bubble injection, [s]

Stop Time

Time to end bubble injection, [s]

Bubble Density

Density of the bubble, [kg/m 3 ]

Bubble Injection Ramp Time

Time over which the bubble injection rate increases from 0 to the resolved feed rate, [s]. Not applicable for Dump or Sink types.

Initial Bubble Size Distribution

This parameter allows the bubble diameter, \(d\), to be randomly sampled from a distribution function at run-time.

Distribution Type

Six initial bubble diameter distributions are available


Single diameter injection.


Diameter randomly chosen from the probability density function defined by min/max, such that:

\[\begin{split}p(d)=\begin{cases}{\frac {1}{\textrm{max}-\textrm{min}}}&{\text{for }}d\in [\textrm{min},\textrm{max}]\\0&{\text{otherwise}}\end{cases},\end{split}\]

where the mean bubble diameter is defined by:


Diameter randomly chosen from the probability density function defined by \(\mu\) and \(\sigma\), such that

\[p(d)=\displaystyle {\frac {1}{\sigma {\sqrt {2\pi }}}}e^{-{\frac {1}{2}}\left({\frac {d-\mu }{\sigma }}\right)^{2}},\]

where the mean bubble size is defined by \(\mu\).


Diameter randomly chosen from the probability density function defined by \(\mu\) and \(\sigma\), such that

\[p(d)=\displaystyle {\frac {1}{d\sigma {\sqrt {2\pi }}}}\ \exp \left(-{\frac {\left(\ln d-\mu \right)^{2}}{2\sigma ^{2}}}\right),\]

where the mean bubble size is defined by:

\[d_m=\displaystyle \exp \left(\mu +{\frac {\sigma ^{2}}{2}}\right).\]

Diameter randomly chosen from the probability density function defined by \(\sigma\), such that:

\[p(d)= \frac{d}{ \sigma ^{2} }e^{-d^{2}/\left(2\sigma^{2}\right)},\]

where the mean bubble size is defined by:

\[d_m=\sigma {\sqrt {\frac {\pi }{2}}}.\]

A comparison between the available probability density functions is shown below. Note that each of these distributions, with the parameters specified in the legend, generates a mean bubble diameter of 0.005 m.


Initial Velocity

Initial Velocity

Initial velocity vector of the bubbles entering the system, [m/s]

Bubble Forces

Bubble trajectories are governed by Newton’s second law, such that:

\[m_i \frac{d \mathbf{v}_i}{dt}=\mathbf{F}_i^f+\mathbf{F}_i^c+\mathbf{F}_i^g\]

where \(m_i\) is the mass of bubble \(i\), \(\mathbf{v}_i\) is the bubble velocity vector, \(\mathbf{F}_i^f\) is the fluid-bubble interaction force, \(\mathbf{F}_i^c\) is a custom user-defined force, and \(\mathbf{F}_i^g\) is the net gravity force on the bubble, which includes the effects of buoyancy.

The fluid-bubble interaction force on bubble \(i\) is given by:


where \(\mathbf{F}_p\) is the pressure gradient force, \(\mathbf{F}_v\) is the virtual mass force, and \(\mathbf{F}_s\) is the Saffman lift force, and \(\mathbf{F}_d\) is the drag force.

The pressure gradient force is given by:

\[\mathbf{F}_p=-V_p \nabla p,\]

where \(V_p\) is the volume of the bubble and \(\nabla p\) is the pressure gradient.

The virtual mass force is given as [1]:

\[\mathbf{F}_m=\left(2.1-\frac{0.132}{0.12+A_c^2} \right) V_p \rho_f \left[ \frac{\left(\dot{\mathbf{u}}-\dot{\mathbf{v}} \right)}{2} \right],\]


\[A_c=\frac{\left(\mathbf{u}-\mathbf{v}\right)^2}{d_p \frac{d \left(\mathbf{u}-\mathbf{v}\right)}{dt}},\]

where \(\rho_f\) is the density of the fluid, \(d_p\) is the bubble diameter, and \(\mathbf{u}\) is the fluid velocity.

The Saffman lift force is given as [2]:

\[\mathbf{F}_s=1.61 d_p^2 \rho_f \sqrt{\frac{ \nu_f }{|\mathbf{\omega}|} } \left[ \left(\mathbf{u}-\mathbf{v}\right) \times \mathbf{\omega_c} \right],\]

where \(\nu_f\) is the kinematic viscosity of the fluid and \(\mathbf{\omega}\) is the fluid vorticity.

The drag force is given by:

\[\mathbf{F}_d=\frac{\pi}{8}C_D d_p^2 \rho_f \left|\mathbf{u}-\mathbf{v} \right| \left(\mathbf{u}-\mathbf{v}\right),\]

where \(C_D\) is the drag coefficient.

Two options are available for defining \(C_D\): “Free Particle” and Packed Bed”

In the Free Particle option, the analytic representation Brown and Lawler [3] describing the drag force on a single bubble in fluid is applied:


In the “Packed Bed” option, the analytic representation of Rong, Dong and, Yu [4] for describing the drag force on a bubble in a packed bed is applied:

\[C_D= \left( 0.63+\frac{4.8}{\sqrt{Re_p}} \right) \epsilon_f^{\left({2-\beta}\right)}\]


\[\beta=2.65 \left( \epsilon_f+1 \right)-(5.3-3.5 \epsilon_f) \epsilon_f^2 e^{\frac{-\left(1.5-\log Re_p^2\right) }{2}},\]

where \(\epsilon_f\) is the fluid volume fraction. In both models, \(Re_p\) is the bubble Reynolds number calculated using the bubble diameter, local fluid viscosity, and the velocity of the bubble relative to the fluid.

Two Way Coupling

Enable/Disable two-way bubble-fluid coupling

Two options are avilable for describing two-way coupling: Density and Third Law.


Under this option, a local body force, \(a_b\), is applied to each fluid voxel in accordance with the local gas-phase volume fraction:

\[a_b=g\alpha \left[ \frac{ \rho_g}{\rho_f}-1\right ],\]

where \(g\) is gravity, \(\alpha\) is the gas-phase volume fraction of the voxel, \(\rho_g\) is the density of the gas phase, and \(\rho_f\) is the density of the liquid phase.

Third Law

The sum of all forces applied to the bubbles in each voxel is applied (in opposite direction) to the fluid.

Bubble Breakup and Coalescence


Two methods are available for modeling bubble break-up and coalescence: (i) Constant Volume Scale and (ii) Dynamic Volume Scale.

Constant Volume Scale

In the discrete approach, every bubble in the gas phase is tacked individually and evolves according to Newton’s second law. Individual bubbles break-up into two individual daughter bubbles following a user-defined break-up expression. Likewise, bubble coalescence occurs when two bubbles collide with an approach Reynolds number, :\(Re_a\), that exceeds a user-set critical value. Break-up and coalescence events explicitly increase and decrease the number of bubbles modeled within the system. As such, the maximum number of bubbles that can be tracked is limited by GPU memory (about 25 million bubbles per GB of RAM).

Dynamic Volume Scale

In the dynamic volume scale approach, the gas phase is modeled as a bubble parcel: a cluster of bubbles defined by a characteristic bubble diameter and a parcel number scale. The characteristic bubble diameter informs the Newtonian trajectory of a parcel; the number scale describes the number of bubbles (with this characteristic diameter) that are represented by the parcel. These two parameters combine to define a total parcel gas surface area and total parcel gas volume. Within this approach, the characteristic diameter adapts dynamically according to a user-defined equilibrium bubble diameter expression. The parcel number scale is then adjusted automatically in a manner than conserves total parcel gas volume during diameter changes. Since the number of parcels modeled within the system does not change due to break-up and coalescence events, this approach has lower memory requirements than the discrete approach.

The difference between the discrete and parcel can be understood as follows: The Contstant Volume Scale appraoch is more accurate, as each bubble is modeled expliclity. With decreasing bubble diameter, however, bubble trajectories become increasingly aligned with the fluid streamlines. The degree of this alignment is characterized by the bubble Stokes number: when the Stokes number is less than 0.1, the difference between the bubble trajectory and the fluid streamline is less than 1%. As bubbles become increasingly aligned with the fluid flow, the necessity of modeling each bubble individually becomes diminished. Rather, as in the parcel model, a representative bubble with an appropriate number scale is sufficient to describe the net-effect of these small gas bubbles on fluid flow and mass transfer. The minimum bubble diameter should therefore be set to a value that corresponds to a low Stokes number. In conventional gas/liquid systems, this condition corresponds to a minimum bubble diameter of about 0.5 mm - 1 mm. Where possible, the Constant Volume Scale should be applied.

Breakup Fraction Expression

For the discrete breakup method, this expression specifies (i) the fluid and bubble properties required initiate a break-up event, and (ii) the volume fraction of the first daughter generated by the breakup process. In order to conserve total volume during the break-up event, the volume fraction of the second daughter generated by the break-up process is automatically defined as the complement of the first volume fraction. In principle, break-up events occur when the local fluid shear forces exceed the local bubble surface tension forces. In practice, a breakup event is assumed to occur if the diameter of a bubble is greater than the equilibrium diameter predicted from turbulence theory. The volume fraction of the daughter bubbles is similarly informed by local fluid shear rates, bubble size, and bubble surface tension. This relationship is defined by the user and can exploit a random number. Note that a volume fraction of 0.5 implies equal-volume daughter bubbles.

Coalescence Reynolds Number

For the discrete breakup method, this expression specifies the critical approach Reynolds number required to initiate a coalescence event. The critical approach Reynolds number is defined as:

\[Re_a=\frac{\mathbf{U}_n d_{ij}}{\nu}\]


\(\mathbf{U}_n\) is the approach velocity of two impacting bubbles, \(d_{ij}\) is the harmonic mean diameter of the bubble pair, and \(\nu\) is the local fluid viscosity. The suggested value for the Coalescence Reynolds Number is 40. [5]

Coalescence User Defined expression

For the discrete breakup method, this user-defined expression specifies the conditions requried for two colliding bubbles to coalese (verus bounce).

Equilibrium Diameter Expression

For the parcel approach this expression defines the local equilibrium bubble diameter as a function of the local fluid properties. The value compute here will define characteristic bubble diameter within the parcel. This expression typically takes the form [6]:

\[D_p=\frac{ \sigma^{\frac{3}{5}} }{ \rho^{ \frac{3}{5} }\epsilon^{ \frac{2}{5} }},\]


\(\sigma\) is the bubble surface tension, \(\rho\) is the density of the fluid surrounding the bubble, and \(\epsilon\) is the local energy dissipation rate.

The bubble diameter equilibration process occurs instantaneously.

Bubble Mass Transfer Coefficient

The mass transfer coefficient, \(k_L\), is a proportionality constant between mass flux and a concentration difference, such that []:

\[\dot{n}=k_L A \Delta C,\]


\(\dot{n}\) is a mass flow rate [mol/s], \(A\) is the mass transfer surface area [m 2 ], and \(\Delta C\) is the driving concentration difference [mol/m 3 ]

The mass transfer coefficient, which has units of m/s, is typically a function of the local fluid properties, such as [6]:

\[k_L=C \left(\epsilon \nu \right)^{a} \textrm{Sc}^{b},\]


\(\epsilon\) is the specific energy dissipation rate [W/kg], \(\nu\) is the fluid kinematic viscosity [m 2 /s ], and \(\textrm{Sc}\) is the fluid Schmidt number.

The local mass transfer coefficient is often multiplied by the local specific bubble surface area (bubble surface area per unit volume) to calculate the local overall mass transfer coefficient, “\(k_La\)”, which has units of [1/s]. Multiplying each local overall mass transfer coefficient by the local concentration difference gives the local mass transfer rate. Integrating this local mass transfer rate over the entire volume gives the total mass transfer rate between the bubbles and the fluid.


The local and total specific surface areas, along with the local and total overall mass transfer coefficients are reported automatically.

Bubble/Scalar Coupling

As discussed in Section Bubble Mass Transfer Coefficient, mass transfer between a bubble and the surrounding fluid is informed by the local fluid mass transfer coefficient, the surface area of the bubble, and the species concentration difference. With the Bubble/Scalar Coupling option, users can model species exchange between the fluid and individual bubbles, along with the associated changes in bubble size and composition (if multiple species are present).

Initial Volume Fraction

This value defines the fraction of the bubble initially occupied by each scalar. If multiple coupled scalar fields are present, the sum of all volume fractions must equal unity.

Molar volume

Volume of space occupied by each mole of gas in the bubble for each species, [m 3/ mol ] This value defines the change in the bubble volume due to mass transfer with the surrounding fluid. A value of zero implies the mass transfer for this species involves no change in the bubble diameter. Note that the standard volume of ideal gas at standard temperature and pressure is 0.0224 m 3/ mol.

Coupling Model

Two methods are available for modeling mass transfer (i) Fick’s Law and (ii) Custom.

Fick’s Law

In this approach, mass transfer is modeled as:

\[\dot{n}=k_L A \left( C_{g,i}-C_{f,i} \right)\]


\(\dot{n}\) is a mass flow rate between the bubble and the surrounding fluid [mol/s], \(k_L\) is the mass transfer coefficient of the fluid surrounding the bubble [m/s], \(A\) is the surface area of the bubble [m 2], \(C_g\) is the saturated concentration of species \(i\) [mol/m 3 ], and \(C_f\) is the local fluid concentration of species \(i\) [mol/m 3 ].

Note that the user must define the saturated concentration of each species.


In this approach, the user defines a custom mass transfer model for each species per local fluid properties and the diameter/composition of each bubble. These models can follow the functional form of Fick’s law, or be extended to model reaction/dissolution processes.


[1] Odar, F. and Hamilton W.S., 1964. Forces on a sphere accelerating in fluid, J. Fluid Mech. 18, 302–314.

[2] Saffman, P.G., 1965. The lift on a small sphere in a slow shear flow, J. Fluid Mech. 22, 385

[3] Brown, P.P. and Lawler, D.F., 2003. Sphere Drag and Settling Velocity Revisited. J. Environ. Chem. Eng. 129, 222-231.

[4] Rong, L.W., Dong, K.J., Yu, A.B., 2013. Lattice-Boltzmann Simulation of fluid flow through packed beds of Uniform Spheres: effect of porosity. Chem. Eng. Sci. 99, 44–58.

[5] Boshenyatov, B., 2012. Laws of Bubble Coalescence and their Modeling. New Developments in Hydrodynamics Research 08, 211-240.

[6] Kewase, Y. and Moo-Young, M., 1990. Mathematical models for design of bioreactors: Applications of: Kolmogoroff’s theory of isotropic turbulence. Chem. Eng. J. 43 B19-B41.

[7] Bird, R.B., Stewart, W.E., Lightfoot, E.N., 2006. Transport Phenomena, 2nd Ed. John Wiley & Sons. New York.