The breakage rate expression, or kernel [400], is expressed as

(14–704)

The birth rate of particles of volume due to breakage is given by

(14–705)

where particles of volume break per unit time, producing particles, of which a fraction represents particles of volume . is the number of child particles produced per parent (for example, two particles for binary breakage).

The death rate of particles of volume due to breakage is given by

(14–706)

The PDF is also known as the particle fragmentation distribution function, or daughter size distribution. Several functional forms of the fragmentation distribution function have been proposed, though the following physical constraints must be met: the normalized number of breaking particles must sum to unity, the masses of the fragments must sum to the original particle mass, and the number of fragments formed has to be correctly represented.

Mathematically, these constraints can be written as follows:

The following is a list of models available in Ansys Fluent to calculate the breakage frequency:

Ansys Fluent provides the following models for calculating the breakage probability density functions (PDF):

The breakage frequency models and the parabolic and generalized PDFs are described in detail in the sections that follow.

The Luo and Lehr models are integrated kernels, encompassing both the breakage frequency and the PDF of breaking particles.

The general form is the integral over the size of eddies hitting the particle with diameter (and volume ). The integral is taken over the dimensionless eddy size . The general form is

(14–710)

where the parameters are as shown in Table 14.4: Luo Model Parameters and Table 14.5: Lehr Model Parameters.

Table 14.4: Luo Model Parameters

	11/3	where	-11/3

For detailed notations, see Luo [400].

Table 14.5: Lehr Model Parameters

	13/3		-2/3

Where is entered through the GUI. See Lehr [353].

The Ghadiri model [201], [460], in contrast to the Luo and Lehr models, is used to model only the breakage frequency of the solid particles. You will have to specify the PDF model to define the daughter distribution.

The breakage frequency is defined as:

(14–711)

where is the particle velocity and is the particle diameter prior to breaking. The breakage constant [s³ m^-11/3] sets the overall scale of how many parent particles is breaking per second. It is case-dependent and is proportional to the particle properties:

(14–712)

where is the particle density, is the elastic modulus of the granule, and is the interface energy. Note that you may need to adjust , for example, by calibrating against a baseline experiment.

The Laakkonen breakage kernel is expressed as the product of a breakage frequency, and a daughter PDF where

(14–713)

where is the liquid phase eddy dissipation, is the surface tension, is the liquid density, is the gas density, is the parent particle diameter and is the liquid viscosity.

The constants , and .

The daughter PDF is given by

(14–714)

where and are the daughter and parent particle volumes, respectively. This model is a useful alternative to the widely used Luo model because it has a simple expression for the daughter PDF and therefore requires significantly less computational effort.

A bubble in a flowing liquid experiences a stress () that deforms the bubble and eventually breaks it once the stress exceeds the restoring effect of the surface tension () as shown in Figure 14.18: Bubble Breakup (adapted from [369]).

Figure 14.18: Bubble Breakup

Based on the work of Martínez-Bazán et al. [416], Liao [369] proposed to define the breakup frequency of a bubble of size in relation to stress and critical stress as follows:

(14–715)

where is the diameter of a child bubble.

The combined breakup approach is used to determine the critical stress based on consideration of energy criterion and force criterion. For binary breakage, where bubble breaks into two bubbles and , can be expressed as:

(14–716)

The energy constraint is determined by the increase of the surface energy during the breakage:

(14–717)

where is the surface energy of the bubble, and the size of bubble is calculated based on the sizes of bubbles and as:

(14–718)

The force constraint is determined by the capillary force of the smallest child bubble:

(14–719)

where is the surface tension coefficient.

The destroying stress that breaks the bubble is caused by the velocity gradient across the bubble’s surface. In a turbulent flow, the velocity variations can result from various mechanisms such as turbulent fluctuation, velocity gradients in a mean flow, turbulent eddies, and interfacial drag [369] as shown in Figure 14.19: Bubble Breakup Mechanisms in a Turbulent Flow (adapted from [369]).

Figure 14.19: Bubble Breakup Mechanisms in a Turbulent Flow

These breakup mechanism stresses are expressed as:

Constants , , , and allow to consider uncertainties in the estimation of destroying stresses caused by different mechanisms. The default (and recommended by Liao [369]) values for the coefficients are: ===1, and = 0.25.

The breakage PDF function contains information on the probability of fragments formed by a breakage event. It provides the number of particles and the possible size distribution from the breakage. The parabolic form of the PDF implemented in Ansys Fluent enables you to define the breakage PDF such that

(14–724)

where and are the daughter and parent particle volumes, respectively. Depending on the value of the shape factor , different behaviors will be observed in the shape of the particle breakage distribution function. For example, if , the particle breakage has a uniform distribution. If , a concave parabola is obtained, meaning that it is more likely to obtain unequally-sized fragments than equally-sized fragments. The opposite of this is true for . Values outside of the range of 0 to 3 are not allowed, because the PDF cannot have a negative value.

Note that the PDF defined in Equation 14–724 is symmetric about .

The generalized form of the PDF implemented in Ansys Fluent enables you to simulate multiple breakage fragments () and to specify the form of the daughter distribution (for example, uniform, equal-size, attrition, power law, parabolic, binary beta). The model itself can be applied to both the discrete method and the QMOM.

Considering the self-similar formulation [539] where the similarity is the ratio of daughter-to-parent size (that is, ), then the generalized PDF is given by

(14–725)

where is the self-similar daughter distribution [141]. The moment of , , is

(14–726)

where

(14–727)

The conditions of number and mass conservation can then be expressed as

(14–728)

(14–729)

The generalized form of [141] can be expressed as

(14–730)

where can be 0 or 1, which represents as consisting of 1 or 2 terms, respectively. For each term, is the weighting factor, is the averaged number of daughter particles, and are the exponents, and is the beta function. The following constraints are imposed on the parameters in Equation 14–730 :

(14–731)

(14–732)

(14–733)

In order to demonstrate how to transform the generalized PDF to represent an appropriate daughter distribution, consider the expressions shown in the tables that follow:

Table 14.6: Daughter Distributions

Type
Equal-size [321]
Attrition [321]
Power Law [674]
Parabolic -a [674]
Austin [34]
Binary Beta -a [259]
Binary Beta -b [420]
Uniform [674]

Table 14.7: Daughter Distributions (cont.)

Type		Constraints
Equal-size [321]
Attrition [321]	2
Power Law [674]
Parabolic -a [674]
Austin [34]
Binary Beta -a [259]	2	N/A
Binary Beta -b [420]	2
Uniform [674]

In Table 14.6: Daughter Distributions, is the Dirac delta function, is a weighting coefficient, and , , , and are user-defined parameters.

The generalized form can represent the daughter distributions in Table 14.6: Daughter Distributions by using the values shown in Table 14.8: Values for Daughter Distributions in General Form.

Table 14.8: Values for Daughter Distributions in General Form

Type									Constraints
Equal-size	1				N/A	N/A	N/A	N/A
Attrition	0.5	2		1	0.5	2	1
Power Law	1			1	N/A	N/A	N/A	N/A
Parabolic	1			2	N/A	N/A	N/A	N/A
Austin				1				1
Binary Beta **	1	2			N/A	N/A	N/A	N/A
Uniform	1		1		N/A	N/A	N/A	N/A

(*)You can approximate by using a very large number, such as 1e20.

(**)Binary Beta -a is a special case of Binary Beta -b when .

The Martinez-Bazan breakage frequency [416] is commonly used to model bubble breakage in stirred tank reactors. The kernel provides a phenomenological expression for the breakage rate of bubbles, which is based on the local values of the turbulence stress and surface tension. The breakage rate is expressed as:

(14–734)

The maximum stable bubble diameter is given by:

(14–735)

In the above equations,

= 0.25 and = 8.2 = constants found experimentally

= parent diameter

= turbulent dissipation

= surface tension

= liquid density

Most daughter distribution PDFs for bubble breakage assume a symmetric U-shaped distribution with the minimum probability of breakage occurring for equi-sized daughter bubbles. Although this kind of a daughter distribution function works well for many applications, an alternative inverted U-shaped distribution has also been proposed in [416], where, by contrast with the U-shaped distribution, the probability of formation of equi-sized daughter bubbles is assumed to be at a maximum. In Ansys Fluent, this is implemented as:

(14–736)

where , , , and is defined in Breakage.

The aggregation kernel [400] is expressed as

(14–737)

where is an aggregation factor for the calibration of the aggregation kernels.

The aggregation kernel has units of m³/s, and is sometimes defined as a product of two quantities:

The birth rate of particles of volume due to aggregation is given by

(14–738)

where particles of volume aggregate with particles of volume to form particles of volume . The factor is included to avoid accounting for each collision event twice.

The death rate of particles of volume due to aggregation is given by

(14–739)

The following is a list of aggregation functions available in Ansys Fluent:

The Luo, free molecular, turbulent, Prince and Blanch, and Liao aggregation functions are described in detail in the sections that follow.

For the Luo model [398], the general aggregation kernel is defined as the rate of particle volume formation as a result of binary collisions of particles with volumes and :

(14–740)

where is the frequency of collision and is the probability that the collision results in coalescence. The frequency is defined as follows:

(14–741)

where is the characteristic velocity of collision of two particles with diameters and .

(14–742)

where

(14–743)

The expression for the probability of aggregation is

(14–744)

where is a constant of order unity, , and are the densities of the primary and secondary phases, respectively, and the Weber number is defined as

(14–745)

Real particles aggregate and break with frequencies (or kernels) characterized by complex dependencies over particle internal coordinates [584]. In particular, very small particles (say up to ) aggregate because of collisions due to Brownian motions. In this case, the frequency of collision is size-dependent and usually the following kernel is implemented:

(14–746)

where is the Boltzmann constant, is the absolute temperature, is the viscosity of the suspending fluid, and are the diameters (or lengths) of particles and . This kernel is also known as the Brownian kernel or the perikinetic kernel.

During mixing processes, mechanical energy is supplied to the fluid. This energy creates turbulence within the fluid. The turbulence creates eddies, which in turn help dissipate the energy. The energy is transferred from the largest eddies to the smallest eddies in which it is dissipated through viscous interactions. The size of the smallest eddies is the Kolmogorov microscale, , which is expressed as a function of the kinematic viscosity and the turbulent energy dissipation rate:

(14–747)

In the turbulent flow field, aggregation can occur by two mechanisms:

For the viscous subrange, particle collisions are influenced by the local shear within the eddy. Based on work by Saffman and Turner [564], the collision rate is expressed as,

(14–748)

where is a pre-factor that takes into account the capture efficiency coefficient of turbulent collision, and is the shear rate:

(14–749)

For the inertial subrange, particles are bigger than the smallest eddy, therefore they are dragged by velocity fluctuations in the flow field. In this case, the aggregation rate is expressed using Abrahamson’s model [6],

(14–750)

where is the mean squared velocity for particle .

The empirical capture efficiency coefficient of turbulent collision describes the hydrodynamic and attractive interaction between colliding particles. Higashitani et al. [245] proposed the following relationship:

(14–751)

where is the ratio between the viscous force and the Van der Waals force,

(14–752)

Where is the Hamaker constant, a function of the particle material, and is the deformation rate,

(14–753)

The aggregation model of Prince and Blanch [535] assumes that the coalescence of two bubbles occurs in three steps:

This aggregation kernel is modeled by a collision rate of two bubbles and a collision efficiency related to the time required for coalescence:

(14–754)

The collision efficiency is modeled by comparing the time required for coalescence and the actual contact time during the collision :

(14–755)

where:

(14–756)

and:

(14–757)

In Equation 14–756 and Equation 14–757, is the density of the liquid, is the surface tension, is the initial film thickness, is the critical film thickness when rupture occurs, is the turbulent eddy dissipation of the liquid, and is the equivalent radius defined as:

(14–758)

where and are the radii of bubbles and , respectively.

The turbulent contributions to collision frequency are modeled as:

(14–759)

where the cross-sectional area of the colliding particles is defined by:

(14–760)

the turbulent velocity is given by:

(14–761)

and is a calibration factor for the turbulent contribution.

The buoyancy contribution to the collision frequency is calculated as:

(14–762)

where:

(14–763)

and is a calibration factor for the buoyancy.

In Equation 14–760, Equation 14–761, and Equation 14–763, and are the diameters of bubbles and , respectively, and is gravity.

The shear contribution to collision frequency ( in Equation 14–754) is currently neglected.

Coalescence occurs where two or more fluid particles collides. This can happen for various reasons such as turbulence, eddy capture, velocity gradient in the mean flow, body forces, and wake entrainment [369] as shown in Figure 14.20: Bubble Collision in a Turbulent Flow (adapted from [369]).

Figure 14.20: Bubble Collision in a Turbulent Flow

Assuming the cumulative contribution of all the coalescence mechanisms, Liao [369] proposed a formulation for the total coalescence frequency of a binary collision between two bubbles of volumes and :

(14–764)

where is the bubble packing limit, and is its maximum. When the local volume fraction of the bubble phase approaches the maximum packing limit , the coalescence rate increases infinitely.

Since bubble-bubble collision can result either in coalescence or bouncing of the colliding bubbles, the coalescence efficiency (or probability) is used. Coalescence frequency is represented as the product of the collision frequency and the efficiency of each collision event :

(14–765)

Similar to the kinetic gas theory, the collision frequency is approximated as the volume swept by the two particles in unit time:

(14–766)

where is the relative velocity between the two bubbles, and is the effective cross-sectional area of collision.

For a random bubble collision caused by turbulent fluctuation, the effective cross-sectional area is:

(14–767)

where and are diameters of bubbles and , respectively.

The relative velocity between bubbles and is estimated by:

(14–768)

where is the turbulent velocity over a distance given by:

(14–769)

where is the turbulent dissipation rate in the liquid phase.

For a turbulence-induced collision, the coalescence frequency between bubble and bubble is obtained as:

(14–770)

where is a constant, and is the coalescence efficiency in the inertial collision regime. Equation 14–770 is valid for bubbles much larger than the Kolmogorov length scale.

When the buoyancy effects are considered, collision can occur only when the faster-moving bubble approaches the slower-moving bubble from bellow in an upward flow. The coalescence frequency is modeled as:

(14–771)

where is a constnat, and and are the terminal velocities of bubble and bubble , respectively. The effective coalescence efficiency is calculated by:

(14–772)

where is the Kolmogorov length scale.

The coalescence efficiency for an inertial collision is expressed as:

(14–773)

where the coefficient has a default value of 5.0, and

(14–774)

where is the liquid density, is the surface tension coefficient, and

(14–775)

The coalescence efficiency of the viscous collision is expressed as:

(14–776)

where is the dynamic viscosity of the liquid, is the Hamaker constant that represents the strength of the van der Waals interactions between macroscopic bodies, and is the characteristic strain rate of the flow in the smallest eddies. For an air-water interface, its value is approximated to 3.7 x 10^-20 J.

In a shear-induced collision, a bubble can only collide with another bubble it follows in the direction of the bulk velocity from the higher-velocity side. The coalescence frequency between bubble and bubble caused by the velocity gradient is expressed as:

(14–777)

where is a constant, is the shear strain rate of the bulk flow, and is defined by Equation 14–772.

For bubbles that are much smaller than the Kolmogorov microscale, the collision is assumed to be controlled by viscous forces. The coalescence frequency is modeled as:

(14–778)

where is a constant. It can be derived from the Kolmogorov scales to be:

(14–779)

When a bubble of a size bigger than a critical diameter moves through a liquid, it generates a wake in the liquid phase behind it. The trailing bubbles within the wake region of the leading bubble can be accelerated and reach the leading bubble. The coalescence efficiency for wake-entrainment is modeled as:

(14–780)

where is a constant, and are the terminal velocities of bubbles and , respectively, and and are the drag coefficients of bubbles and , respectively. The blending function for bubble is defined as:

(14–781)

where the critical diameter is computed as:

(14–782)

where is the gas density.

The constants , , , , and account for uncertainty in the estimation of relative velocity under different conditions. Liao [369] recommends setting these coefficients to 1.0, which are the default settings in Ansys Fluent.

In Ansys Fluent, the PBE is written in terms of volume fraction of particle size :

(14–783)

where is the density of the secondary phase and is the volume fraction of particle size , defined as

(14–784)

where

(14–785)

and is the volume of the particle size . In Ansys Fluent, a fraction of , called , is introduced as the solution variable. This fraction is defined as

(14–786)

where is the total volume fraction of the secondary phase.

The nucleation rate appears in the discretized equation for the volume fraction of the smallest size . The notation signifies that this particular term, in this case , appears in Equation 14–783 only in the case of the smallest particle size.

The growth rate in, Equation 14–783 [256], is discretized as follows:

(14–787)

The volume coordinate is discretized as [256] where and is referred to as the “ratio factor”.

The particle birth and death rates are defined as follows:

(14–788)

(14–789)

(14–790)

(14–791)

where and

(14–792)

is the particle volume resulting from the aggregation of particles and , and is defined as

(14–793)

where

(14–794)

If is greater than or equal to the largest particle size , then the contribution to class is

(14–795)

The default breakage formulation for the discrete method in Ansys Fluent is based on the Hagesather method [224]. In this method, the breakage sources are distributed to the respective size bins, preserving mass and number density. For the case when the ratio between successive bin sizes can be expressed as where , the source in bin , () can be expressed as

(14–796)

Here

(14–797)

A more mathematically rigorous formulation is given by Ramkrishna [326], where the breakage rate is expressed as

(14–798)

where

(14–799)

The Ramkrishna formulation can be slow due to the large number of integration points required. However, for simple forms of , the integrations can be performed relatively easily. The Hagesather formulation requires fewer integration points and the difference in accuracy with the Ramkrishna formulation can be corrected by a suitable choice of bin sizes.

The SMM approach is based on taking moments of the PBE with respect to the internal coordinate (in this case, the particle size ).

Defining the ^th moment as

(14–800)

and assuming constant particle growth, its transport equation can be written as

(14–801)

where

(14–802)

(14–803)

(14–804)

(14–805)

is the specified number of moments and is the nucleation rate. The growth term is defined as

(14–806)

and for constant growth is represented as

(14–807)

Equation 14–802 can be derived by using

and reversing the order of integration. From these moments, the parameters describing the gross properties of particle population can be derived as

(14–808)

(14–809)

(14–810)

(14–811)

(14–812)

These properties are related to the total number, length, area, and volume of solid particles per unit volume of mixture suspension. The Sauter mean diameter, , is usually used as the mean particle size.

To close Equation 14–801, the quantities represented in Equation 14–802 – Equation 14–805 need to be expressed in terms of the moments being solved. To do this, one approach is to assume size-independent kernels for breakage and aggregation, in addition to other simplifications such as the Taylor series expansion of the term . Alternatively, a profile of the PSD could be assumed so that Equation 14–802 – Equation 14–805 can be integrated and expressed in terms of the moments being solved.

In Ansys Fluent, an exact closure is implemented by restricting the application of the SMM to cases with size-independent growth and a constant aggregation kernel.

The quadrature approximation is based on determining a sequence of polynomials orthogonal to (that is, the particle size distribution). If the abscissas of the quadrature approximation are the nodes of the polynomial of order , then the quadrature approximation

(14–813)

is exact if is a polynomial of order or smaller [140]. In all other cases, the closer is to a polynomial, the more accurate the approximation.

A direct way to calculate the quadrature approximation is by means of its definition through the moments:

(14–814)

The quadrature approximation of order is defined by its weights and abscissas and can be calculated by its first moments by writing the recursive relationship for the polynomials in terms of the moments . Once this relationship is written in matrix form, it is easy to show that the roots of the polynomials correspond to the eigenvalues of the Jacobi matrix [534]. This procedure is known as the moment inversion algorithm, which uses either the Product-Difference algorithm [213] or the Wheeler algorithm [705]. Once the weights and abscissas are known, the source terms due to coalescence and breakage can be calculated and therefore the transport equations for the moments can be solved.

Applying Equation 14–813 and Equation 14–814, the birth and death terms in Equation 14–801 can be rewritten as

(14–815)

(14–816)

(14–817)

(14–818)

Theoretically, there is no limitation on the expression of breakage and aggregation kernels when using QMOM.

The nucleation rate is defined in the same way as for the SMM. The growth rate for QMOM is defined by Equation 14–806 and represented as

(14–819)

to allow for a size-dependent growth rate.

The DQMOM equations, describing a poly-dispersed particle system undergoing aggregation, breakage, and growth can be written as follows (the details of the DQMOM formulation can be found in [169]):

(14–820)

(14–821)

where and are the VOF and the effective length of the particle phase, respectively. is the number of particles per unit volume and is the growth rate at Quadrature point , while and can be computed through a linear system resulting from the moment transformation of the particle number density transport equation using Quadrature points. The linear system can be written in matrix form as

(14–822)

Where the coefficient matrix is defined by

(14–823)

(14–824)

The vector of unknowns is defined by

(14–825)

The right hand side of Equation 14–822 is the known source terms involving aggregation and breakage phenomena only. The growth term is accounted for directly in Equation 14–820 and Equation 14–821. At present, nucleation is not considered.

(14–826)

The source term for the moment is defined as

(14–827)

When the abscissas of the Quadrature points are distinct, the matrix is well defined and a unique solution of Equation 14–822 can be obtained. Otherwise, the matrix is not full rank and cannot be inverted to find a unique solution for . The method adopted by Ansys Fluent to overcome this problem is to employ a perturbation technique. For example, for the current three Quadrature points system, the perturbation technique will add a small value to the abscissas to make sure the matrix A is full rank. It is important to note that the perturbation technique is only used for the definition of matrix and no modifications are made to the source term vector of Equation 14–827. Therefore, both the weights and overall source terms resulting from aggregation and breakage are not affected by the perturbation method. The simulation tests have found that the perturbation method can stabilize the solutions of Equation 14–822 and reduce the physically unrealistically large source terms for the two phases whose abscissa are too close in value. However, the technique has little effect on the phase whose abscissa is distinct from the other two.