# Particles¶

## 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 (i) massless tracers that follow the fluid streamlines, (ii) inertial particles with a density/volume, (iii) discrete element method (DEM) particles that interact with the fluid and other DEM particles, and (iv) resolved particles that are large compared to the lattice spacing.

Particles enter the system along the children geometry of the particle parent with a user-defined injection velocity. Particle translational (and rotational, when relevant) 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.

Note

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

## General¶

- Injection Geometry
This option describes where particles 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 particle addition/removal rate.

**Sink**– Remove all particles that enter the geometry**Source**– Maintain a minimum user-defined particle population inside the children geometry over a period**Dump**– Dump a user-defined number of particles into the children geometry at a given time**Feed**– 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- Specified Feed Rate
Only presented for Volume Feed injection options. This is the number of particles 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 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 Option
Only presented for Source, Feed, and Volume Feed injection options. This option defines when to stop adding particles

**End of Simulation**– Inject particles for the entire simulation**Specified**– Stop injection at a specifed user time, [s]- Particle Type
Four different particle types are available

**Tracer**– Particles have no inertia and act as fluid tracers. No particle-particle interactions are considered.**Inertial**– Spherical particles with have a density and diameter that translate according to Newton’s second law. No particle-particle interactions are considered. Fluid/particle force calculated from empirical correlations. Particles are assume to be small relative to lattice spacing.**Discrete Element**– Superquadratic particles with a density and shape that translate/rotate according to Newton’s second law. Particle-particle interactions are considered. Fluid/particle force calculated from empirical correlations. Particles are assume to be small relative to lattice spacing.**Resolved**– Superquadratic particles with a density and shape that translate/rotate according to Newton’s second law. Particle-particle interactions are considered. Fluid/particle force calculated directly from the solver. Particles are assume to be large relative to lattice spacing.- Particle Density
Only presented for Inertial, Discrete Element, and Resolved particle options. This value represents the density of the particle, [kg/m

^{3}]

## 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** – 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 particle diameter of 0.005 m.

## Particle Forces¶

This section is only presented for *Inertial* and *Discrete Element* particle types. Particle trajectories are governed by Newton’s second law, such that:

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:

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

The virtual mass force is given as [1]:

with

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]:

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

The drag force is given by:

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:

with

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¶

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

**Disabled** – Particles do not inform fluid dynamics.

**Third Law** – The sum of all forces applied to the solid particles in each voxel is applied (in opposite direction) to the fluid.

**Density** – A local body force, \(a_b\), is applied to each fluid voxel in accordance with the local solid-phase volume fraction

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

The **Density** approach is a simpler two-way coupling algorithm that only considers the effects of buoyancy/gravity when calculating the local fluid body force. That is, the body force on the fluid voxel is only a function of the particle 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 particles with particle diameter \(d_b\), particle density \(\rho_b\), particle area \(A_b\), and particle 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 particles are gravity/buoyancy \(F_b\) and drag \(F_d\), the body force \(F_T\) on the fluid in the voxel is calculated from:

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

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

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

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

where \(\epsilon\) is the solid-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 particles in agitated tanks. For larger particles moving through quiescent tanks, however, the complete form of Newton’s third law should be considered.

## Static Body Interactions¶

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

For *Tracer* and *Inertial* particles, options include *Bounce* and *Stick*. The *Bounce* option implies that particles bounce off the solid body family. The *Stick* option implies that particles stick to the solid body family.

For Discrete Element particles, options include the *Bounce Simple*, *Bounce Custom*, and *Stick*. Under the “Bounce Simple* option, particle-wall interactions are assumed to have the same interaction parameters as the particle-particle interactions (e.g. the same Young modulus, friction, etc.). Under the *Bounce Custom* option, particle-wall interaction can be assigned custom custom interaction parameters.

For *Resolved* particles, only the *Bounce* option is available.

## Solid Interactions¶

This section is only presented for *Discrete Element* and *Resolved* particle types. These parameters defined how particle interact with other particles within the particle family.

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:

and

where \(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:

and

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

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

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:

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

and

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

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

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

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

subject to:

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

The rolling friction torque is given by:

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

- Track Rotation
This options activates tracking of particle rotation and angular velocity.

## Particle Breakup and Agglomeration¶

Two representations are available for modeling changes in particle size due to break-up and coalescence: (i) Discrete and (ii) Local Equilibrium

- Discrete
In this approach, every particle in the solid phase is tacked individually and evolves according to Newton’s second law. Individual particle break-up into two individual daughter particles following a user-defined break-up expression. Likewise, particle coalescence occurs when two colliding particles satisfy a user-defined coalescence model. 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).

- Breakup Fraction Expression
For the explicit approach, this expression specifies (i) the fluid and particle 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 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 particles 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.

- Coalescence Model
For the explicit approach, this is the algorithm that will determine if two colliding particles bounce off each other or coalesce into one particles. This expression can be a user-defined function or a evoke the colliding particle 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 \(Re_a=\frac{\mathbf{U}_n d_{ij}}{\nu}\) , where \(\mathbf{U}_n\) is the approach velocity of two colliding 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]**Coalescence User Defined expression**– This user-defined expression specifies the fluid and particle pair conditions required for two colliding particles to coalesce (versus bounce).- Coalescence Sampling Interval
Number of time steps to wait between particle collision checks.

- Local Equilibrium
In this 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 particle surface area and total parcel particle volume. Within this approach, the characteristic diameter adapts dynamically according to a user-defined equilibrium particle diameter expression. The number scale of the parcel is then adjusted automatically in a manner than conserves total parcel particle 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 particle diameter as a function of the local fluid properties. The value computed here will define the characteristic particle 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, \(\sigma\) is the effective 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.

Note

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

## Particle Mass Transfer Coefficient¶

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

where \(\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]:

where \(\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)\), where \(\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.**Custom**– 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.

## Advanced¶

- 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. 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 captures 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 particle count.- Initial Velocity
Initial velocity vector of the particles entering the system, [m/s]

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

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

## Define discrete superquadrics (SQ) size distribution¶

Edit the discrete superquadrics (SQ) size distribution brings up the editor dialog. Edit the table to define the SQ size and shape parameters. The probabilities must sum up to 1.0.

**Probability**– The probability that this particle will be created**a**– The radius in the x direction. Length units shown in the current unit system of the model.**b**– The radius in the y direction. Length units shown in the current unit system of the model.**c**– The radius in the z direction. Length units shown in the current unit system of the model.**n1**– Shape parameter 1**n2**– Shape parameter 2

## References¶

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