Basic Concepts

Particles are Lagrangian spheres that can access continuous positions across the fluid lattice. They are used to model two-phase flow, slurry transport processes, dissolution processes, and suspensions. Particles can either be massless tracers that follow the fluid streamlines or as particles with mass/inertia.

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

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

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

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


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



Name of the particle set


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

Injection Geometry

Two options are available

Inject at surface

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

Inject in volume

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

Injection Option

Six different injection options are available


Remove all particles that enter the geometry


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


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


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

Volume Feed

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

Injection Number Scale

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

For large-scale systems with high volume feed rates, which may have billions of particles, the computational requirements necessary to monitor every particle may exceed available memory. To make such simulations tractable, users may specify an Injection Number Scale that reduces the number of explicitly tracked particles by this factor. The simulated volume of the particles, 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 particles that will be fed into the system, [#/s]

The resolved feed rate describes the rate at which particles 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 particle volume. Otherwise, the resolved feed rate will be smaller by a factor equal to the Injection Number Scale. Note that particles are injected uniformly across the entire child(ren) geometry surface area (or volume, if injecting in the child volume).

Start Time

Time to begin particle injection, [s]

Stop Time

Time to end particle injection, [s]

Particle Density

Density of the particle, [kg/m 3 ]

Particle Injection Ramp Time

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

Track Rotation

For inertial particles, this options activates tracking of particle rotation and angular velocity.

Initial Particle Size Distribution

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

Distribution Type

Six initial particle 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 particle 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 particle 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 particle 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 particle 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 particle diameter of 0.005 m.


Initial Velocity

Initial Velocity

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

Particle Forces

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

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

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

The fluid-particle interaction force on particle \(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 particle 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 particle 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 particle 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 particle 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 particle Reynolds number calculated using the particle diameter, local fluid viscosity, and the velocity of the particle relative to the fluid.

Two Way Coupling

Enable/Disable two-way particle-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 solid-phase volume fraction:

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

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

Third Law

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

Particle-Particle Interactions

Interactions between two particles \(i\) and \(j\) are described by the Hertz contact model [8]:


where \(f_{n,ij}\) and \(f_{t,ij}\) are the normal and tangential forces, defined by:





\(k_{n,ij}\) and \(k_{t,ij}\) are the normal and tangential stiffnesses, \(\gamma_{n,ij}\) and \(\gamma_{t,ij}\) are the normal and tangential damping coefficients, \(\delta_{n,ij}\) and \(\delta_{t,ij}\) are the normal and tangential overlaps, and \(\dot{\delta_{n,ij}}\) and \(\dot{\delta_{t,ij}}\) are the time derivatives.

The stiffness are related to physical properties of the particles and the overlap distances, such that:

\[k_{n,ij}=\frac{4}{3} Y^*_{ij} \sqrt{R^*_{ij} \delta_{n,ij} },\]


\[k_{t,ij}= 8 G^*_{ij} \sqrt{R^*_{ij} \delta_{n,ij} },\]

with the equivalent Young’s modulus, \(Y^*_{ij}\), defined by:

\[\frac{1}{Y^*_{ij}}= \frac{1-\nu_i^2}{Y_i}+\frac{1-\nu_j^2}{Y_j},\]

and the equivalent shear modulus, \(G^*_{ij}\), defined by:

\[\frac{1}{G^*_{ij}}= 2\frac{\left(2+\nu_i\right)\left(1-\nu_i\right)}{Y_i}+2\frac{\left(2+\nu_j\right)\left(1-\nu_j\right)}{Y_j},\]

where \(Y\), \(G\), and \(\nu\) are the Young’s modulus, Shear modulus, and Poisson ratio of particle \(i\) and the equivalent radius, \(R^*_{ij}\), is defined by:

\[\frac{1}{R^*_{ij}}= \frac{1}{R_i}+\frac{1}{R_j}.\]

The damping parameters are likewise related to physical properties of the particles and the overlap distances:

\[\gamma_{n,ij}=-2\left(\sqrt{\frac{5}{6}}\right)\left(\frac{\ln{e_r}}{\sqrt{\ln^2{e_r}+\pi^2}}\right)\sqrt{\frac{2}{3} k _{n,ij} m^*_{n,ij} },\]


\[\gamma_{t,ij}=-2\left(\sqrt{\frac{5}{6}}\right)\left(\frac{\ln{e_r}}{\sqrt{\ln^2{e_r}+\pi^2}}\right)\sqrt{k_{t,ij} m^*_{n,ij} },\]

where \(m^*_{ij}\), in addition to making for a terrific company name, is the equivalent mass defined by:

\[\frac{1}{m^*_{ij}}= \frac{1}{m_i}+\frac{1}{m_j},\]

where, \(Y\), \(G\), and \(\nu\) are the Youngs modulus, Shear modulus, and Poisson ratio of particle :math:i`.

When particle rotation is tracked, the time-evolution of the angular velocity of particle \(i\), \(\omega_i\), follows Newtons second law:

\[I_i \frac{d \mathbf{\omega}_i}{dt}=\displaystyle \sum_j \left(\mathbf{M}_{t,ij}+\mathbf{M}_{r,ij} \right),\]

where \(j\) are the particles interacting with particle \(i\), and the tangential torque is given by:

\[\mathbf{M}_{t,ij}=\mathbf{r}_{i} \times f_{ct,ij},\]

subject to:

\[\left| f_{ct,ij} \right| \leq \mu_{s,ij} \left| f_{cn,ij} \right|,\]

where \(\mathbf{r}_i\) is the distance to the contact point.

The rolling friction torque is given by:

\[\mathbf{M}_{r,ij}=-\mu_{r,ij} \left| f_{cn,ij} \right| \frac{\mathbf{\omega}_{ij}}{\left|\mathbf{\omega}_{ij}\right|}R^*_{ij},\]

where \(\mu_{s,ij}\) and \(\mu_{r,ij}\) are the sliding and rolling friction coefficients between the particle pair.

Particle Breakup and Agglomeration


Two methods are available for modeling particle break-up and coalescence: (i) discrete and (ii) parcel

Discrete breakup

In the discrete approach, every particle in the solid phase is tacked individually and evolves according to Newton’s second law. Individual particles break-up into two individual daughter particles following a user-defined break-up expression. Likewise, particle coalescence occurs when two particles 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 particles modeled within the system. As such, the maximum number of particles that can be tracked is limited by GPU memory (about 25 million particles per GB of RAM).

Parcel breakup

In the parcel approach, the solid phase is modeled as a particle parcel: a cluster of particles defined by a characteristic particle diameter and a parcel number scale. The characteristic particle diameter informs the Newtonian trajectory of a parcel; the number scale describes the number of particles (with this characteristic diameter) that are represented by the parcel. These two parameters combine to define a total parcel solid surface area and total parcel solid volume. Within this approach, the characteristic diameter adapts dynamically according to a user-defined equilibrium particle diameter expression. The parcel number scale is then adjusted automatically in a manner than conserves total parcel solid 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.

Breakup Fraction Expression

For the discrete breakup method, this expression specifies (i) the fluid and particle 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 particle surface tension forces. In practice, a breakup event is assumed to occur if the diameter of a particle is greater than the equilibrium diameter predicted from turbulence theory. The volume fraction of the daughter particle is similarly informed by local fluid shear rates, particle size, and particle 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 particles.

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

Agglomeration Reynolds Number

For the discrete breakup method, this expression specifies the critical approach Reynolds number required to initiate a agglomeration 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 particles, \(d_{ij}\) is the harmonic mean diameter of the particle pair, and \(\nu\) is the local fluid viscosity. The suggested value for the Coalescence Reynolds Number is 40. [5]

Equilibrium Diameter Expression

For the parcel and hybrid approaches, this expression defines the local equilibrium particle diameter as a function of the local fluid properties. The value compute here will define characteristic particle 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 particle surface tension, \(\rho\) is the density of the fluid surrounding the particle, and \(\epsilon\) is the local energy dissipation rate.

The particle diameter equilibration process occurs instantaneously.

Particle 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 particle surface area (particle 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 particles and the fluid.

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

Particle/Scalar Coupling

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

Initial Volume Fraction

This value defines the fraction of the particle 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 solid in the particle for each species, [m 3/ mol ] This value defines the change in the particle 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 particle diameter.

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 particle and the surrounding fluid [mol/s], \(k_L\) is the mass transfer coefficient of the fluid surrounding the particle [m/s], \(A\) is the surface area of the particle [m 2], \(C_g\) is the saturated concentration of species \(i\) [mol/m:sup: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 particle. These models can follow the functional form of Fick’s law, or be extended to model reaction/dissolution processes.


Compute Particle 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 particles across the fluid domain to that of random particle distribution. In this sense, the comparison can be used as an estimate on particle 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 particles per unit volume. As the particles homogenize through the system, the average nearest neighbor distance should approach this value. A average nearest neighbor distance below this value implies particle 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 particles and \(N\) is the total number of particles in the system.


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

[8] Blais, B., Bertrand, F., 2017. CFD-DEM investigation of viscous solid–liquid mixing: Impact of particle properties and mixer characteristics. Chem. Eng. Res. Des. 118, 270–285.

[9] Chandrasekhar, S., 1943. Stochastic Problems in Physics and Astronomy. Rev. Mod. Phys. 15, 1–89.