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 can enter the system along the children geometry of the bubble parent or through an inlet boundary condition. Bubble positions and velocities evolve according to Newton’s Second Law using a Verlet algorithm advanced at the simulation timestep. As discussed in Section Fluid Interaction 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 image files are printed at the same frequency as slices and volumes.

General Properties

Injection Geometry

This option describes where bubbles will enter the system

  • Child surface - Bubbles will be injected along the surface of the children geometry within the bubble family

  • Child volume - Bubbles will be injected within the volume enclosed by the children geometry with the bubble family

  • Boundary Condition - Bubbles will be injected at a user-selected inlet.

Injection Option

This option characterizes the bubble addition/removal rate.

  • Sink - Remove all bubbles that enter the geometry.

  • Source - Maintain a minimum user-defined bubble population inside the children geometry or inlet over a period.

  • Dump - Inject a user-defined number of bubbles into the children geometry or inlet as a function of time.

  • Feed - Continuously feed bubbles into the children geometry or inlet at a user-defined number flow rate over a period.

  • Volume Feed - Continuously feed bubbles into the children geometry or inlet at a user-defined volume flow rate as a function of time. If the Bubble Pressure option is on, the bubble density is defined by the reference pressure and temperature (See Section Advanced).

  • Specified Feed Rate - Only presented for Volume Feed injection options. This is the number of bubbles that will be fed into the system, [#/s] When the Injection Number Scale is set to 1 (See Section Advanced), 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 Option

Only presented for Source, Feed, and Volume Feed injection options. This option defines when to stop adding bubbles

  • End of Simulation - Inject bubbles for the entire simulation

  • Specified - Stop injection at a specified user time, [s]

Bubble Density

Density of the bubbles, [kg/m 3 ]. If the Bubble Pressure option is on, the bubble density is defined by the reference pressure and temperature (See Section Advanced).

Initial Bubble Size Distribution

Distribution Type

Six initial bubble diameter probability distribution functions are available. At runtime, the bubble diameter will be randomly sampled from the given distribution. If the Bubble Pressure option is on, the bubble diameter is defined by the reference pressure and temperature (See Section Advanced).

Single – Single diameter injection.

Uniform – Diameter randomly chosen from the probability density function defined by min/max, such that: \(p(d)=\begin{cases}{\frac {1}{\textrm{max}-\textrm{min}}}&{\text{for }}d\in [\textrm{min},\textrm{max}]\\0&{\text{otherwise}}\end{cases}\), where the mean particle diameter is defined by: \(d_m=\tfrac{1}{2}(\textrm{max}+\textrm{min})\).

Normal – 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 particle size is defined by \(\mu\).

Log-normal – 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 particle size is defined by, \(d_m=\displaystyle \exp \left(\mu +{\frac {\sigma ^{2}}{2}}\right).\)

Rayleigh – 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 particle size is defined by: \(d_m=\sigma {\sqrt {\frac {\pi }{2}}}.\)

Discrete – Diameter randomly chosen from a discrete probability mass function defined by the user. Note that probabilities associated with each possible diameter must be positive and sum to unity

A comparison between the five continuous 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.


Fluid Interaction

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.

Pressure Gradient

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.

Virtual Mass

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.

Saffman Lift

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.

Drag Force

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.

Drag Force Option

Three options are available for defining \(C_D\): Disable, Free Particle, and Packed Bed.

  • Disabled - No drag force on the bubbles.

  • Free Particle - The analytic representation Brown and Lawler [3] describing the drag force on a single bubble in fluid is applied.

  • Packed Bed - 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 the Free Particle and Packed Bed 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

This setting affects how the bubbles inform fluid flow. Three options are available (i) Disabled, (ii) Third Law and (iii) Density.


Bubbles do not inform fluid dynamics.


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 fluid voxel hosting the bubble, \(\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.

The Density approach is a simpler two-way coupling algorithm that only considers the effects of buoyancy when calculating the local fluid body force. That is, the body force on the fluid voxel is only a function of the bubble volume fraction within the voxel. Although the Third Law option presents a more rigorous application of the conservation of momentum, small time steps may be required to keep the simulation stable if large body forces are realized in the fluid.

To assess the suitability of the “Density” approach, consider a single fluid voxel containing \(N\) identical bubbles with bubble diameter \(d_b\), bubble density \(\rho_b\), bubble area \(A_b\), and bubble volume \(V_b\). Take the continuous phase fluid density to be \(\rho_f\) and the fluid voxel volume to be \(V_l\).

If the primary forces on the bubbles are buoyancy \(F_b\) and drag \(F_d\), the body force \(F_T\) on the fluid in the voxel is calculated from:

\[F_T=F_b+F_d=\sum_i^N [g V_b \left(ρ_b-ρ_f\right)+\frac{1}{2} ρ_f A_b C_d V_s^2 ] ,\]

where \(C_d\) is the drag coefficient on the bubbles, \(g\) is acceleration due to gravity, and \(V_s\) is the bubble slip velocity, and the sum is over all bubbles in the voxel. Since bubbles are assumed to be identical, this expression can be re-cast into:

\[F_T=g N_b V_b (ρ_b-ρ_f )[1+\lambda]\]

where \(N_b\) is the number of bubbles in the voxel and

\[\lambda=\frac{ρ_f A_b C_d V_s^2}{2g V_b \left(ρ_b-ρ_f\right) }.\]

Recognize that, for spherical gas bubbles in turbulent liquids with order-of-magnitude larger densities, this expression can be simplified to:

\[\lambda \approx \frac{V_s^2}{gd}.\]

For fluids, we are interested primary in the volumetric body force, such that

\[f_T=\frac{F_T}{V_b} =g \epsilon \left(ρ_b-ρ_f\right)\left[1+\lambda\right]\]

where \(\epsilon\) is the gas-phase volume fraction within the voxel that varies between 0 and 1.

The \(\lambda\) term represents the ratio between drag force and buoyancy force. When :math:`lambda’ is much smaller than 1, the body force on the voxel due to buoyancy is much greater than the body force on the fluid voxel due to fluid drag. As, such the volume fraction approach to two-way coupling is sufficient. This condition is typically realized for millimeter sized bubbles in agitated tanks. For larger bubbles moving through quiescent tanks, however, the complete form of Newton’s third law should be considered.

Custom Fluid Bubble Force

User defined expression for an additional external force applied the bubbles. This expression can be a function of bubble properties and local fluid properties.


Bubble position and velocities are updated in tandem with the fluid flow using the simulation timestep.

Static Body Interactions

This parameter specifies how each bubble set interacts with each solid body family.

Options include Bounce and Stick. The Bounce option implies that bubbles bounce off the solid body family. The Stick option implies that bubbles stick to the solid body family.

Bubble Breakup and Coalescence

Two algorithms are available for modeling changes in bubble size due to break-up and coalescence: (i) Explicit and (ii) Parcel


In this 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 colliding bubbles satisfy a user-defined coalescence model. 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).

Breakup Fraction Expression

For the explicit approach, this expression specifies (i) the fluid and bubble properties required to 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 Model

For the explicit approach, this is the algorithm that will determine if two colliding bubbles bounce off each other or coalesce into one bubble. This expression can be a user-defined function or a evoke the colliding bubble Reynolds number.

Critical Reynolds

This value defines the critical approach Reynolds number required to initiate a coalescence event. The critical approach Reynolds number is defined as: .. math:

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

where \(\mathbf{U}_n\) is the approach velocity of two colliding 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

This user-defined expression specifies the fluid and bubble pair conditions required for two colliding bubbles to coalesce (versus bounce).

Coalescence Sampling Interval

Number of time steps to wait between bubble collision checks.

Local Equilibrium

In this 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 number scale of the parcel is then adjusted automatically in a manner than conserves total parcel gas volume when the 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.

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 computed here will define the characteristic bubble diameter within the parcel. This expression typically takes the form [6]:

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

where \(C\) is a constant between 0.1-1.0, \(\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.


For appropriate sampling, the number of individually tacked parcels should be within an order-of-magnitude of the number of fluid voxels.

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 local specific energy dissipation rate [W/kg], \(\nu\) is the local fluid kinematic viscosity [m 2 /s ], and \(\textrm{Sc}\) is the local 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 between the bubble and it’s surrounding fluid gives the mass transfer rate for the bubble. 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 be equal to one.

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 volume of ideal gas at standard temperature and pressure is 0.0224 m 3/ mol. If the Bubble Pressure option is on, the bubble density is defined by the reference pressure and temperature (See Section Advanced).

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 the 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 more complex reaction/dissolution processes.


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. Note that if break-up and coalescence are active and the Injection Number Scale is greater than 1, the local equilibrium representation should be applied as not all coalescence will be captured explicitly.

Resolved Feed Rate

Shown for reference as part of the Volume Feed injection option. This is calculated from the actual feed rate divided by the injection scale. This value is useful when estimating the system bubble count.

Compute Bubble Distribution Data

Compute the average nearest neighbor separation distance and nearest neighbor separation distance probability distribution function. These metrics can be used to compare the spatial arrangement of bubbles across the fluid domain to that of random bubble distribution. In this sense, the comparison can be used as an estimate on bubble homogenization Note that, for a random distribution of hard spheres, the average distance between neighbors, \(a\), can be estimated as:

\[a=\left({\frac {3}{4\pi n}}\right)^{1/3}\]

where \(n\) is the number of bubbles per unit volume. As the bubbles homogenize through the system, the average nearest neighbor distance should approach this value. A average nearest neighbor distance below this value implies bubble clumping.

The nearest neighbor separation distance probability distribution function, \(P_n\), is then defined from [9]:

\[P_n(r)=\frac{3}{ a}\left(\frac{r}{a}\right)^2 \left(1 - \left(\frac{r}{a}\right)^3 \frac{1}{N} \right)^{N - 1}\]

where \(r\) is the distance between bubbles and \(N\) is the total number of bubbles in the system.

Deletion Diameter

Bubbles with diameters below this critical value will be deleted from the simulation. Any residual concentration will be dumped into the local fluid voxel.

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.

Bubble Pressure

If active, the effects of local hydrostatic pressure on bubble size are considered. This functionality is only applicable to free surface simulations. Per the ideal gas law, the bubble density is directly related to pressure while the molar volume is inversely related to pressure.

The reference pressure, which is used as a basis for the volume flow rate, initial bubble size, and molar volume, is presented. By default, the reference pressure is taken to be the initial pressure of the headspace gas.


[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.