3.3.1. Agglomerator

The new method of solving the agglomeration PBE was implemented as a new model of the Dyssol simulation system according to Equation (29). To simplify the model, it was assumed that the mass within the apparatus remains constant and there is no classification by size at the output. Therefore, the mass flow of material leaving the apparatus . *mout*(*t*) equals to the mass flow entering it . *min*(*t*), and the distribution of the material at the output *Y*(*t*) equals to the distribution in the holdup *<sup>H</sup>*(*t*).

Since the applicability of the proposed method does not depend on the chosen agglomeration kernel, to describe the agglomeration rate in all case studies, a frequently used physically relevant kernel based on a Brownian motion [55,56] was applied in the following from:

$$\begin{aligned} \boldsymbol{\beta}(\mathbf{x}, \mathbf{x}') &= \beta\_0 \boldsymbol{\beta}^\* \Big( \mathbf{x}^{\frac{1}{3}} + \mathbf{x}'^{\frac{1}{3}} \Big) \mathbf{x}^{-\frac{1}{3}} + \mathbf{x}'^{-\frac{1}{3}} \Big), \\ \boldsymbol{\beta}^\* &= \begin{cases} 1, & \mathbf{x} + \mathbf{x}' \le \mathbf{x}\_{\mathbf{z}, \min} \\ 1 - \frac{(\mathbf{x} + \mathbf{x}') - \mathbf{x}\_{\mathbf{z}, \min}}{\mathbf{x}\_{\mathbf{z}, \max} - \mathbf{x}\_{\mathbf{z}, \min}}, & \mathbf{x}\_{\mathbf{z}, \min} < \mathbf{x} + \mathbf{x}' < \mathbf{x}\_{\mathbf{z}, \max} \\ 0, & \mathbf{x} + \mathbf{x}' \ge \mathbf{x}\_{\mathbf{z}, \max} \end{cases} \end{aligned} \tag{57}$$

Here β(*<sup>x</sup>*, *x* ) denotes the agglomeration rate between particles *x* and *x* ; β0 is the size-independent agglomeration rate constant, dependent on operating conditions; β∗ is the size- dependent agglomeration rate constant; and *xa*,*min* and *xa*,*max* are minimum and maximum sizes of resulting agglomerates.

The agglomerator unit was implemented according to Equations (28) and (29) to calculate and apply the transformation matrix for the particles size. To consider secondary distributions, it was extended to perform aggregation with averaging of the properties distributed over the second dimension, which is calculated as follows.

Let particles of two sizes *xsrc*1 and *xsrc*2 aggregate into a particle of size *xdst*. Then *Tsrc*1,*dst* and *Tsrc*2,*dst* are the corresponding entries of the transformation matrix, calculated by Equation (28). The summarized mass fraction of particles from size-class *i* is distributed over *L*2 classes of the secondary dimension, so that

$$c\_i = \sum\_{j=1}^{L\_2} c\_{i,j\prime} \tag{58}$$

where *ci*,*j* is a mass fraction of granular material with size *xi* and value of the secondary dimension falling into class *j*. If particles from two classes (*src*1, *j*1) and (*src*2, *j*2) aggregate, the result lands into some class (*dst*, *j*). The size-class *dst* is directly determined from the agglomeration mechanism and obtained from the transformation matrix. The mass fraction of each *j*–th class of the secondary distribution is calculated as

$$c\_{\rm dst,j} = \sum\_{j=1}^{L\_2} \sum\_{j2=1}^{L\_2} \frac{(c\_{\rm src2,j,1} T\_{\rm src1,dst}) \cdot \left( c\_{\rm src2,j2} \cdot T\_{\rm src2,dst} \right)}{\sum\_{k=1}^{S\_2} c\_{\rm src2,k} \cdot T\_{\rm src2,dst}} \tag{59}$$

whereas the class index *j* for equidistant grid can be obtained as

$$j = \frac{j1 \cdot c\_{\rm src1,j1} \cdot T\_{\rm src1,dst} + j2 \cdot c\_{\rm src2,j2} \cdot T\_{\rm src2,dst}}{c\_{\rm src1,j1} \cdot T\_{\rm src1,dst} + c\_{\rm src2,j2} \cdot T\_{\rm src2,dst}}.\tag{60}$$

3.3.2. Mill

The mill model, developed in Dyssol, implements Equations (47) and (48) to calculate transformation laws and Equations (52) and (55) to apply them. To calculate Equations (41) and (43), the adaptive Simpson's method of numerical integration [57] was applied. The criterion used to determine when to terminate interval subdivision was chosen according to Lyness [58], and its desired tolerance was set to 10−15.

For simplification, the holdup mass within the unit stays constant and no classification by size occurs at the outlet.

As can be seen from Section 2.3, the proposed approach does not depend on the chosen selection or breakage functions. Therefore, for testing purposes, the following physically relevant models were taken from literature. To describe the selection rate in all case studies, the selection function based on the model proposed by King [59] was applied:

$$S(\mathbf{x}) = s\_0 \cdot \begin{cases} 0, & \mathbf{x} \le \mathbf{x}\_{b,\text{min}} \\ 1 - \left(\frac{\mathbf{x}\_{b,\text{max}} - \mathbf{x}\_{b,\text{min}}}{\mathbf{x}\_{b,\text{max}} - \mathbf{x}\_{b,\text{min}}}\right)^n, & \mathbf{x}\_{b,\text{min}} < \mathbf{x} < \mathbf{x}\_{b,\text{max}} \\ 1, & \mathbf{x} \ge \mathbf{x}\_{b,\text{max}} \end{cases} \tag{61}$$

where *<sup>S</sup>*(*x*) is the mass fraction of particles of size *x* that are crushed; *s*0 introduces a size-independent selection rate factor; *xb*,*min* and *xb*,*max* are critical sizes of particles; and *n* is a power law exponent of the King's selection function.

The breakage function was implemented according to Vogel et al. [60] in a number-based distributed form: 

$$B(\mathbf{x}, y) = \begin{cases} 0.5 \frac{q}{y} \Big(\frac{\mathbf{x}}{y}\Big)^{q-2} \Big(1 + \tanh\left(\frac{y-y'}{y'}\right)\Big), & y \ge \mathbf{x} \\\ 0, & y < \mathbf{x} \end{cases}.\tag{62}$$

Here *<sup>B</sup>*(*<sup>x</sup>*, *y*) denotes the mass fraction of particles with the size of *y*, which after breakage ge<sup>t</sup> the size less or equal to *x*; *y* is the minimum fragment size that can be achieved after milling; and *q* is a power law exponent of the Vogel breakage function.
