Morphology-dependent random binary fragmentation of in silico fractal-like agglomerates

Fragmentation is important in nature and industry. Polymer degradation [1], fracture of solids and volcanic eruptions [2], network topology and resilience [3], droplet [4] and agglomerate [5] breakup in turbulent flow are just a few examples where fragmentation determines the outcome of a process.

The fragmentation of linear chains has been extensively studied in polymer science since bond degradation or depolymerization breaks up the polymer. In most cases random scission [1,6,7], whereby all bonds break with equal probability, is considered. In colloids, the size distribution of colloidal aggregates in a fluid flow, and thence the suspension rheology, is greatly influenced by hydrodynamic shear-induced stresses that may break or restructure the particles [8].

In aerosol science, which plays a fundamental role in, e.g., atmospheric sciences (cloud formation), air pollution (ultimate fate of atmospheric particulate pollutants), and industrial production of pharmaceuticals, the particle size distribution is one of the most important parameters that specifies the chemical and physical properties of the aerosol.

The dynamics of the size distribution, as determined by agglomeration, whereby the colliding particles retain their identity (or coagulation whereby colliding particles coalesce) and fragmentation, is usually described by mean-field rate equations. These population balance equations are generalizations of the Smoluchowski coagulation equation [9]. The kinetics of irreversible agglomeration has been extensively analysed, often with emphasis on scaling properties of the size distribution or the characteristic cluster size [10]. As clusters grow their number decreases continuously, their characteristic size increases, and the size distribution shrinks. In contrast, if fragmentation occurs, a likely event as the ever-growing clusters are more likely to fragment, a steady-state distribution may be established in the long-time limit.

Our primary interest in this work lies in the characterization of the size distribution of fragments (and the associated fragmentation kernel) upon random scission of fractal-like agglomerates, and its dependence on the morphology of the initial agglomerate, as described by its fractal dimension (for fixed fractal prefactor). We consider numerically the linear, binary fragmentation of synthetic, loop-less, fractal-like agglomerates upon random link removal. Agglomerates that satisfy a scaling law over a finite range (fractal-like) are generated in silico. They are composed of equal-size, point-touching spherical monomers.

Consider the irreversible fragmentation equation [11], for the discrete particle size distribution n_i (t), the number of agglomerates composed of i monomers (of size i) per unit volume at time t,

$$\begin{equation} \frac{\mathrm{d} n_i (t)}{\mathrm{d}t} = \sum_{j = i + 1}^{i_{\textrm{max}}} a_j p_{ij} n_j(t) - a_i n_i(t), \end{equation} \tag{ 1 }$$

where a_i is the (herein taken to be time-independent) fragmentation rate of a particle of size i, and p_ij is the distribution of fragments of size i resulting from the break-up of a particle of initial size j. The maximum particle size $i_{\textrm{max}}$ is usually taken to be infinity.

The fragment size distribution must obey a number of conservation laws [11,12]. The first constraint gives the expected number of fragments upon single-bond removal. For binary fragmentation it reduces to

$$\begin{equation} \sum_{i=1}^{j-1} p_{ij} = 2, \quad \quad \forall j. \end{equation} \tag{ 2 }$$

Since each particle breaks into two fragments, the fragment size distribution must be symmetric, $p_{ij} = p_{(j-i)j}$. The second conservation law is mass conservation: the sum of the number of monomers in all the fragments must equal the number of monomers in the initial agglomerate,

$$\begin{equation} \sum_{i=1}^{j-1} i \, p_{ij} = j, \quad \quad \forall j. \end{equation} \tag{ 3 }$$

In appendix A of the Supplementary Material Supplementarymaterial.pdf (SM) we show that for binary fragmentation eq. (2) implies eq. (3), rendering only one conservation law independent. In addition, eq. (2) suggests that $g_{ij} = p_{ij}/2$ defines a (discrete) probability distribution. Then, eq. (3) shows that the average fragment size is $j/2$ (as expected).

The fragmentation rate a_j is estimated by assuming that all bonds break with equal probability. Let the lifetime of each bond be characterized by a probability distribution, e.g., a Poisson distribution with average bond lifetime (characteristic time of fragmentation) τ. The number of bonds in a j-mer that would break in time $t = \tau$ is the number of monomers times the probability of a bond breaking. Thus, the fragmentation frequency is proportional to the initial size of the agglomerate, $a_j = j/\tau$.

It is occasionally convenient to approximate the discrete particle size distribution by a continuous distribution. This approximation is very accurate for $j \gg 1$. We also introduce a vector of morphological descriptors $\vec{t}{s}$ to account for the morphology of the aggregate [13]. Accordingly, the discrete size distribution n_i(t) is replaced by the continuous distribution $n(x,t; \vec{t}{s})$ where $n(x,t; \vec{t}{s}) \text{d} x$ is the number of particles per unit volume in the size range x and $x + \text{d} x$ and at morphological state $\vec{t}{s}$. For the cases considered herein, the additional morphological descriptor is the fractal dimension $\vec{t}{s} = d_f$. The discrete variables i and j are replaced by the continuous variables x (fragment size) and y (initial agglomerate size), respectively. The continuous version of eq. (1) becomes

$$\begin{eqnarray} \frac{\mathrm{d} n (x,t;d_f)}{\mathrm{d}t} &= \int_{x}^{{x}_{\textrm{max}}} \text{d} y a(y) p (x,y; d_f) n(y,t; d_f)\nonumber\\ &\quad {} - a(x) \, n(x,t; d_f), \end{eqnarray} \tag{ 4 }$$

where p(x, y; d_f) is the fragmentation density function to obtain a x-mer from the fragmentation of a y-mer. It is related to its discrete counterpart via $p(x,y; d_f) \text{d} x = p_{ij}$. The fragmentation frequency is $a(x) = x/\tau$. The fragmentation density function must satisfy the conservation laws (see, also, appendix A of the SM),

$$\begin{equation} \int_0^y \text{d} x p(x,y; d_f) = 2, \qquad \int_0^y \text{d} x x p(x,y; d_f) = y \end{equation} \tag{ 5 }$$

with the symmetry condition $p(x,y; d_f) = p(y-x,y; d_f)$. As in the discrete case, only one conservation law is independent. The function p(x, y;d_f) is defined for $1 \le x \le y-1$. Henceforth, we will present results for the continuous fragment size distribution. The discrete case is summarized in appendix C of the SM.

We generated in silico independent, fractal-like agglomerates of specified structure (fixed fractal prefactor $k_f = 1.3$ and varying dimension $d_f = 1.6, 1.8, 2.1$) and mass (number of monomers j). We used the hierarchical tunable, off-lattice, cluster-cluster agglomeration algorithm, initially proposed by Thouy and Jullien [14] and later modified by Filippov et al. [15], to create mostly loop-less, agglomerates composed of equal-sized, spherical, point-touching monomers. We only considered loop-less agglomerates. The algorithm does not aim to reproduce a physical agglomeration mechanism in that the motion of the colliding clusters is not explicitly simulated. Instead, it is based on geometric (including steric) considerations. The advantage of such a geometric algorithm is that the desired cluster morphology (d_f, k_f) is an input (whereas in explicit cluster simulations it is an output). The synthetic agglomerates satisfy by construction exactly the fractal-like scaling law ${j} = k_f (R_g/R_1)^{d_f}$, where R_g is the radius of gyration of the cluster [16] and R₁ is the monomer radius. As summarized in appendix B, table I, of the SM the number of agglomerates varied from 400000 to 25000. The number of monomers per cluster varied from approximately 80 to 800.

According to the cluster generation algorithm (appendix B of the SM) clusters are created via monomer-monomer contact. If the resulting cluster does not have loops, the addition of a cluster onto an existing cluster generates only one bond. Thus, the total number of bonds in a loop-less j-cluster is j –1; the sum of bonds connecting all monomers in a cluster (over both pairs $(i,j-i)$ and $(j-i,i)$) is $2(j-1)$. If f_ij is the number of bonds that upon breaking would divide the initial cluster into two fragments of sizes i and j –i, then their sum is $\sum_{i=1}^{j-1} f_{ij} = 2 (j-1)$, as discussed extensively in ref. [17]. The corresponding fragment distribution becomes $p_{ij} = f_{ij}/(j-1)$. Equivalently, these clusters have a (mean) coordination number of $c_j = 2(j-1)/j$ [ [16,18].

The fragmentation algorithm is based on mapping the initial agglomerate onto a symmetric adjacency matrix. Each monomer configuration (i.e., cluster) has a unique representation as a graph via the adjacency matrix. A symmetric adjacency matrix A_ij identifies completely a non-directed graph. The adjacency matrix is used to determine the connected components of a structure, thence the fragments upon removal of a bond. A similar approach was used by Isella and Drossinos [19] to identify clusters generated via monomer Langevin dynamics.

The adjacency matrix is constructed from all monomer-monomer Euclidean distances: monomer-monomer links (bonds) are represented by one and their absence by a zero. Two monomers (i, j) are considered bonded if their center-of-mass distance D_ij is smaller than the monomer diameter 2 R₁ plus a threshold distance $D_{\textrm{thr}}$, i.e., $D_{ij} < 2 R_1 + D_{\textrm{thr}}$. We took $D_{\textrm{thr}} = 10^{-3} \times 2R_1$. A pictorial representation of a fragmentation event (referred to as "neck breakage") via adjacency matrices is shown in fig. 1.

Zoom In Zoom Out Reset image size

The fragmentation algorithm consists of randomly choosing a non-zero element of the adjacency matrix according to a discrete uniform probability distribution, $\xi \epsilon U(0,N_{\textrm{link}})$, where $N_{\textrm{link}}$ is the number of links in the agglomerate. The selected element is replaced by a zero, thereby eliminating the monomer-monomer bond. A new adjacency matrix is formed, and its connected components are determined. If the initial agglomerate fragments, the sizes of the fragments are stored. If the aggregate does not fragment, it is removed from the list of initial clusters. Only one link is removed per agglomerate, no multiple fragmentation events are considered. Then a new cluster is chosen, and the same procedure is repeated till the list of clusters has been exhausted. We found that fewer than 0.05% of the clusters did not fragment. By removing them we ensure that we consider only loop-less agglomerates.

Figure 2 shows two simulated fragmentation events. In the left panel the resulting fragments are dissimilar, on the right they are of equal size. We found that the fragments have a slightly smaller fractal dimension, consequently their structure is more open, whereas they have a slightly higher fractal prefactor: they tend to be more locally compact than the initial agglomerate [16]. The structural parameters of the two fragments were not identical: they showed a slight asymmetry, which disappears in the fragment distribution.

Zoom In Zoom Out Reset image size

Figure 3 shows representative histograms of the fragment size distribution upon single-bond random fragmentation of small (left, j = 96) and large (right, j = 770) Diffusion Limited Cluster Aggregation (DLCA) agglomerates $(d_f = 1.8$, $k_f = 1.3)$. It is apparent that this fragmentation algorithm leads to largely dissimilar fragment sizes: the resulting fragment size distribution is U-shaped. Kalay and Ben-Naim [20] also found a U-shaped fragment distribution in their study of the fragmentation of random trees by removing a single node (a slightly different fragmentation process). A predominance of non-equal–sized fragments upon bond removal in random clusters on a lattice was also found in ref. [3]. The experimental results of Kusters et al. [21], who found that ultrasonic fragmentation of agglomerated particles in liquids leads to erosion, namely only a limited number of monomers break up from the agglomerate, also suggest a U-shaped fragment size distribution.

Zoom In Zoom Out Reset image size

We fitted the empirically determined fragment distribution to a symmetric beta distribution because it reproduces a U-shaped distribution for certain values of its parameter. Of course, our choice is not unique: other distributions could reproduce our data, for example, eq. (11). However, the beta-distribution fit is more parsimonious: it depends on one parameter only (for fixed number of monomers in the agglomerates and fixed aggregate morphology (d_f, k_f)). Moreover, a symmetric beta distribution was introduced ad hoc previously to model the fragment size probability density function [22], see also ref. [23]. The functional form of the fragment density function p(x, y), expressed in terms of the relative number of monomers $z= x/y$, was chosen to be

$$\begin{equation} p(x,y; d_f) = \frac{2}{y} b[z; \beta(d_f, y)], \end{equation} \tag{ 6a }$$

with

$$\begin{equation} b[z; \beta(d_f, y)] = \frac{[z(1-z)]^{\beta(d_f,y)}} {\int_{z_{\text{min}}}^{z_{\text{max}}} \, \text{d} z \, [z(1-z)]^{\beta(d_f,y)}} \end{equation} \tag{ 6b }$$

for

$$\begin{equation} \frac{1}{y} \equiv z_{\text{min}} \leq z \leq 1-\frac{1}{y} \equiv z_{\text{max}} \end{equation} \tag{ 6c }$$

By construction the continuous probability density function b(z) is symmetric, $b(z) = b(1-z)$, and it satisfies the conservation laws

$$\begin{eqnarray} &&\int_{z_{\text{min}}}^{z_{\text{max}}} \text{d} z \, b[z; \beta(d_f,y)] = 1, \end{eqnarray} \tag{ 7a }$$

$$\begin{eqnarray} &&\int_{z_{\text{min}}}^{z_{\text{max}}} \text{d} z \, z \, b[z; \beta(d_f,y)] = \frac{1}{2}. \end{eqnarray} \tag{ 7b }$$

$[i^{-1}+{(j-i)}^{-1}]^\alpha$

$[z(1-z)]^{-\alpha}$

For a given initial cluster size y, and a fractal dimension, the fragment distribution may be obtained by minimizing the distance, in the sense of the l²-norm, between the numerical results (the histograms shown in fig. 3 appropriately normalized) and the predictions of eqs. (3). The minimization procedure would minimize the distance between $b(z) \text{d} z$ and $p_{ij}/2$ for all $i = 1, \ldots, j-1$ to ensure that $b(z) \text{d} z = p_{ij}/2$. For j –1 fragments $\text{d} z = 1/j$, and the optimization condition becomes $b(z) = j p_{ij}/2$. This procedure leads to a global optimization, i.e., an optimization over all fragment sizes.

Instead, we opted for a local optimization. It is more important to reproduce the small-fragment behaviour $({z \to 1/y}$ or $z \to 1-1/y)$, the size of most fragments, than the complete distribution. Moreover, since the number of equal-size fragments is relatively small, the quality of the simulation data deteriorates as $z \to 1/2$. We determined the probability of occurrence of a monomer fragment for all the clusters we generated for the three fractal dimensions. We found that

$$\begin{equation} g_{1y} = \frac{1}{2} p_{1y} = -0.1442 + 0.1518 d_f, \qquad \forall y, \end{equation} \tag{ 8 }$$

provides an excellent fit of our simulation data. Equation (8) suggests that p_1y depends only on the fractal dimension and not on the size of the initial agglomerate. The exponent β was obtained by requiring that the probability of monomer occurrence

$$\begin{equation} g_{1y} = \frac{1}{y} b\left[\frac{1}{y}; \beta(d_f,y)\right] = \frac{1}{y} \, \frac{[\frac{1}{y}(1-\frac{1}{y})]^{\beta(d_f,y)}} {\int_{z_{\text{min}}}^{z_{\text{max}}} \, \text{d} z \, [z(1-z)]^{\beta(d_f,y)}}, \end{equation} \tag{ 9 }$$

equal the empirically determined monomer occurrence probability. We used eq. (8) to generate monomer occurrence probabilities for $j =10, 11, \ldots, 1010$, and we fitted the resulting beta-distribution exponent to

$$\begin{equation} \beta(d_f,y) = a(d_f) + b(d_f) \, y^{c(d_f)} \end{equation} \tag{ 10a }$$

to obtain the coefficients

$$\begin{eqnarray} a(d_f) & = -0.42 d_f - 0.44, \quad b(d_f) = -1.28 d_f + 12.29, \nonumber\\ c(d_f) & = -0.37 d_f - 0.31, \quad {\text{\textrm{for}} \quad d_f \epsilon [1.6, 2.1]}. \end{eqnarray} \tag{ 10b }$$

We show the exponent $\beta(d_f,y)$ in fig. 4. It tends to a constant asymptotic value as $y \to \infty$. In that limit the fragment density function p(x, y; d_f) becomes approximately a homogeneous function under the scaling $(x,y) \rightarrow (\lambda x, \lambda y)$ since $\lim_{y \to \infty} \beta(d_f,y) = \beta(d_f)$. The kernel does not become exactly homogeneous because some y-dependence remains in the integration limits. The y-dependent limits are required because the asymptotic exponent is ${\beta(d_f, \infty) < -1}$ (for the studied fractal dimensions). The beta distribution is properly defined only for ${\beta > -1}$, the reason why our empirical fit eq. (6) is defined for $z \, \epsilon \, (0,1)$. This is not a limitation since there is always a finite monomer size [23]. Homogeneous fragmentation kernels have been used extensively in the past, see, e.g., refs. [11,12,25,26]. We found that the coefficients presented in eq. (10b) depend slightly on whether a discrete or continuous distribution was modelled. The differences are small, but noticeable; see fig. 1 in appendix C of the SM.

Zoom In Zoom Out Reset image size

Predictions of our empirical fit are compared to numerical simulations of single-bond, random fragmentation of DLCA agglomerates in fig. 5. Simulation results (filled square red symbols) were obtained by normalizing (to unity) the histogram of fragment sizes with j –1 bins. The beta-distribution points, converted from the continuous density function to a discrete mass function, are denoted by filled blue pentagrams (note the logarithmic y-scale of the main figure; inset: double logarithmic plot).

Zoom In Zoom Out Reset image size

We also compare the numerical distributions to an expression proposed by Odriozola et al. [17], who analysed the fragment distribution of DLCA clusters generated by a slightly different algorithm. They fitted the fragment size distribution, specifically, f_ij which is the number of bonds that upon breaking would divide a cluster into fragments of sizes i and j –i, to a four-parameter function,

$$\begin{equation} {f}_{ij} = p_1 [ i^{p_2} + {(j-i)}^{p_2}] [ i^{-p_3} + {(j - i)}^{-p_3}] [ i {(j-i)}]^{p_4}, \end{equation} \tag{ 11 }$$

with $p_1 = 0.4391$, $p_2 = 1.006$, $p_3 = 1.007$, and $p_4 = -0.1363$. The product of the second term times the third term is reminiscent of the continuum Brownian agglomeration kernel for two fractal-like clusters of dimension d_f. Note that $p_2 \sim p_3 \sim 1$, suggesting that eq. (11) approximates a symmetric beta distribution with a fixed exponent, independent of j, and for a fixed fractal dimension $(d_f = 1.8)$.

As argued earlier, eq. (11) may by normalized to unity by dividing it by the total number of bonds to obtain

$$\begin{equation} {g}_{ij} = \frac{{f}_{ij}}{2 (j-1)}. \end{equation} \tag{ 12 }$$

The predictions of eq. (12) are denoted by filled, green circles in figs. 5,6, and 7. We note reasonable agreement of simulation data with theoretical predictions. It should be stressed that the normalization condition $\sum_{i=1}^{j-1} {f}_{ij}$ is approximately equal to $2({j}-1)$. In fact, the authors of ref. [17] used this condition to estimate the accuracy of their fit. This implies that eq. (12) does not ensure mass conservation, rendering its use in population balance equations problematic. This is in sharp contrast to eqs. (6) that by construction satisfy exactly mass conservation.

The effect of morphology is studied in fig. 6. The two empirical fits are compared graphically to numerically determined fragment distributions arising from the break-up of 192-monomer clusters for the three fractal dimensions. As the fractal dimension increases the probability of obtaining small fragments increases, the fragment distribution sharpens. As before the agreement is very good, even though eq. (11) was derived for relatively small DLCA clusters (j < 100).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

We also present the probability of occurrence of a monomer fragment as a function of initial agglomerate size parametrized by the fractal dimension (of the initial agglomerate) in fig. 7. Since the exponent β was determined from the monomer probability, our proposed fit reproduces very well the simulation data, and it shows a dependence on the fractal dimension. As the fractal dimension decreases, the probability of monomer occurrence decreases, In fact, we expect that as the fractal dimension decreases the fragment distribution would flatten out: the straight-chain limit $(d_f = 1)$ is a constant probability distribution independent of the fragment size, but dependent on y, $b(z;d_f=1) =1, p(x,y;1) = 2/y$.

We also studied the effect of morphology through a universal fragment size probability density function (a parameter-free distribution) that depends only on the initial cluster fractal dimension. We constructed it by requiring that it satisfy exactly the conservation laws, eqs. (7), and that it have the correct limit for straight chains. The universal fragment size probability distribution is

$$\begin{eqnarray} b_{\textrm{uni}}(z; d_f) &= \Big \{ d_f z(1-z) [ z^{1/d_f} + (1-z)^{1/d_f}]\nonumber\\ &\quad {}\times[ z^{-1/d_f} + (1-z)^{-1/d_f}]\Big \}^{-1}. \end{eqnarray} \tag{ 13 }$$

It is easy to show that the conservation law is satisfied: $\int_0^1 \text{d} z\, b_{\textrm{uni}}(z;d_f ) = 1$ as calculated analytically by Mathematica [27]. In the straight-chain limit $b_{\textrm{uni}}(z, d_f = 1) = 1$ ($p(x,y; d_f = 1) = 2/y$). Moreover, since eq. (13) is a function of z only and independent of y, under the scaling of sizes, $b_{\textrm{uni}}(z;d_f)$ remains invariant. Its homogeneity index is zero, leading to a homogeneous fragmentation kernel. The predictions of eq. (13) are compared to those of the beta distribution in fig. 8, top panel (for the asymptotic beta-distribution exponent $\beta(d_f,1000)$). It is apparent that the parameter-free distribution does not reproduce the simulation data for $d_f \neq 1$ as accurately as the beta distribution. However, it does reproduce the straight-chain limit, bottom panel. That panel shows that the effect of the fractal dimension becomes significant only for fractal dimensions very close to the singular limit $d_f \to 1$. The transition from a pronounced U-shaped distribution to the limiting case of a constant (flat) distribution of a straight line becomes noticeable for d_f < 1.1 in that at $d_f=1.1$ the probability of a monomer fragment occurring is still predicted to be almost twice that of a monomer fragment of a straight chain.

Zoom In Zoom Out Reset image size

We studied numerically the fragment size distribution upon random scission of in silico fractal-like agglomerates as a function of the initial size of the agglomerate and its structural properties, fractal dimension $(d_f = 1.6, 1.8, 2.1)$ and fixed prefactor $(k_f = 1.3)$. The fragment size distribution is an essential ingredient of the fragmentation kernel. Earlier studies, for example refs. [10,12,24,26], used fragmentation kernels in population balance equations whose analytical form was dictated either by physical arguments or by properties of related homogeneous agglomeration kernels, instead of originating from simulation data or experimental measurements. We followed the latter approach to obtain an analytical expression for the fragment size distribution that reproduces rather accurately our numerical simulations.

Our numerical simulations showed that, for our fragmentation algorithm, the fragment size distributions are U-shaped, namely most clusters fragment into largely dissimilar fragments. We showed that a symmetric beta distribution reproduces rather accurately the empirical fragment size distributions. We found that the beta-distribution exponent depends on the initial agglomerate morphology (via the fractal dimension for fixed prefactor) and the number of monomers y in the initial agglomerate, tending to a d_f-dependent asymptotic limit for large y.

We, also, derived a universal (parameter-free) fragment size distribution, dependent on the fractal dimension of the initial cluster, by requiring that it satisfy the two fragmentation conservation laws and that it reproduce the straight-chain limit.

The views expressed are purely those of the authors and may not in any circumstances be regarded as stating an official position of the European Commission. ADM was partially supported by the Horizon 2020 EU Framework Programme through the SUREAL-23 project (Grant Agreement 724136).