**1. Introduction**

Bulk materials typically consist of individual non-uniform particles having different parameters, which vary in a certain range. These parameters are referred to as distributed and include, for example, the diameters of particles, their densities, or porosities. Although these parameters can physically be fairly accurately described by continuous distribution functions, such a representation is not always convenient for numerical analysis and modelling. Therefore, in practice, the continuously distributed material parameter space is discretized into shorter intervals, usually called classes. Each class is assigned a quantity of material, whose parameters fall into this particular interval. In this case, the distribution of the parameter of the entire material in this class is narrowed down from a continuous distribution function to a single value, for example, the average value of this class. Thus, it is assumed that all the particles inside have the same value of the parameter. Therefore, a simulation system can operate using a finite number of parameter values, which (if the number of classes is sufficient) accurately describe the original distribution.

Despite the fact that most models in solids processing technology consider particle size distribution as one of the main material parameters, many of them additionally incorporate various other attributes. For example, the shape and orientation of primary particles are important in crystallization processes [1]; the size and yield strength of the colliding granules can be considered as parameters for wet granulation models [2]; the distributions of granules by porosity and saturation are important for breakage rate in some apparatuses [3]; and the moisture content of particles plays a significant role in dryers [4].

Such a diversity of models and their parameters makes the simulation process much more challenging when trying to combine their different types in a single flowsheet. The main problem here is that each unit is usually being developed to consider only a limited number of distributed parameters which are important for its model, neglecting others. Moreover, the usual approach to develop models involves a direct calculation of the distributed parameters at the outlet of a unit. As a result, if the material is additionally described by some secondary distributed parameters, information about them and their dependencies may be lost during the simulation.

For example, consider a flowsheet where a screen and a dryer are connected in series, and the solid phase is distributed according to particle size and moisture content. In this case, the screen unit, usually performing separation of particles only according to size, will fail to determine their moisture properly, which is important for the dryer.

To handle this issue, Pogodda [5] introduced an approach using transformation matrices and implemented it into the SolidSim simulation environment. Later, Dosta [6] and Skorych et al. [7] extended this technique for the dynamic flowsheet simulation using the Euler method. Transformation matrices in this approach describe the laws of mass transfer between all classes of the distribution. Their application, instead of explicit calculation of output flows, enables the usage of more material parameters in a proper way and implicitly preserves information about all secondary distributions and their dependencies at any time point in the unit.

#### *1.1. Flowsheet Simulation of Solid Phase Processes*

Nowadays, flowsheet simulation software is commonly used to aid with design, planning, optimization, and testing in chemical engineering [8]. Usually it can help with the minimization of operation costs, troubleshooting, increasing throughput and product quality, and generally provides a better understanding of the process. However, despite the fact that flowsheet simulation of fluid processes has already been state of the art for many decades [9], generally applicable systems for modelling granular materials began to appear only recently and are now being actively developed. This arises from the complexity of description and handling of the solid phase. In contrast to fluids, it requires treatment of distributed parameters, such as particle size, internal porosity or moisture content. Moreover, these distributions can be interrelated, which requires even more advanced methods to treat them properly and preserve their interconnection during the simulation.

Currently, flowsheet modelling systems that can work with the solid phase are being actively developed. Among them may for example be mentioned:


Active development of new dynamic models and usage of flowsheet simulations to study various solid phase processes, such as continuous granulation [15,16], chemical looping combustion [17], crystallization [1], milling processes [18], separation [19], tablet manufacturing [20,21], and agglomeration [22] are of increasing interest in this area.

#### *1.2. Use of Population Balance Equations for Particulate Processes*

Population balance equations (PBE) appear in various branches of engineering where particles continuously change their properties, like in pharmaceutical industries and in the processing of minerals, food, paints, and ink. For example, they are used to describe time-dependent processes in crystallizers, fluidized beds, liquid–liquid, gas–liquid contactors, or polymer reactors. In many particulate processes, one-dimensional PBEs with particle size<sup>1</sup> as the property coordinate are used. In this article, the size is referred to as the volume of particles, unless explicitly stated otherwise. These equations describe a time-dependent change of particle size distribution (PSD) in a specific volume, which is caused by various processes, such as agglomeration, breakage, nucleation, growth, or shrinkage. This contribution considers only agglomeration and breakage processes. The time evolution of the particle size distribution *<sup>u</sup>*(*<sup>x</sup>*, *t*) ≥ 0 in these processes, at time *t* ≥ 0, for particles of size *x* ≥ 0 can be represented by a special type of nonlinear integro-di fferential equation, namely a population balance equation [23,24]

$$\begin{split} \frac{\partial u(\mathbf{x},t)}{\partial t} &= \frac{1}{2} \int\_{0}^{\infty} \beta(\mathbf{x}', \mathbf{x} - \mathbf{x}') u(\mathbf{x}', t) u(\mathbf{x} - \mathbf{x}', t) d\mathbf{x} \nu - \int\_{0}^{\infty} \beta(\mathbf{x}, \mathbf{x}') u(\mathbf{x}, t) u(\mathbf{x}', t) d\mathbf{x}' \\ &+ \int\_{\mathbf{x}}^{\infty} b(\mathbf{x}, \mathbf{y}) S(y) u(\mathbf{y}, t) dy - S(\mathbf{x}) u(\mathbf{x}, t) + \dot{Q}\_{\text{in}}(t) - \dot{Q}\_{\text{out}}(t). \end{split} \tag{1}$$

Here, the first and the third terms of the right hand side of Equation (1) represent the birth events (appearance of new particles) due to agglomeration and breakage, respectively. The death events (disappearance of existing particles) owing to agglomeration and breakage process are described by the second and the fourth terms, respectively. Terms . *Qin*(*t*) and . *Qout*(*t*) denote the input and the output in a continuous process, accordingly. The function β(*x* , *x*) is the agglomeration kernel, which represents the rate at which a particle of size *x* agglomerates with a particle of size *x*. The breakage function *b*(*<sup>x</sup>*, *y*) in the third term represents the fraction of fragments of size *x*, which appear when a particle of size *y* breaks. The selection function *<sup>S</sup>*(*x*) describes the rate at which a particle of size *x* is selected to break.

In general, the breakage function is described by the following dependencies:

$$\begin{cases} \, \, b(\mathbf{x}, y) = 0, & \forall \mathbf{x} \ge y \\ \, \, \int\_0^y b(\mathbf{x}, y) d\mathbf{x} = \, \, \theta(y), & \forall y > 0 \end{cases} \tag{2}$$

and

$$\int\_{0}^{y} x b(x, y) dx = y, \quad \forall y > 0. \tag{3}$$

Equation (2) designates that the breakage function *b*(*<sup>x</sup>*, *y*) will be equal to zero when the size of the resulting particle *x* is greater than the size of the initial particle *y*. When the initial particle of size *y* disintegrates to a particle of size *x* during the breakage process, the total number of particles can be represented by the function <sup>ϑ</sup>(*y*) and, obviously, <sup>ϑ</sup>(*y*) ≥ 1. Moreover, when a particle of size *y* splits into daughter particles, it generally satisfies the local mass conservation law, which is given by Equation (3).

Various numerical methods have been applied to solve Equation (1). Examples for the most commonly used methods are: the method of moments [25–28], the fixed pivot technique [29,30], the Monte Carlo simulation method [31–34], the finite element method [35–37], approaches based on a separable approximation of the kernel function with the fast Fourier transformation [38–40], the cell average technique [41,42], and the finite volume method [43,44]. In this work, the finite volume method [43,44] is applied for the spatial discretization, and the second-order Runge–Kutta method

is used for the time variable to ge<sup>t</sup> the full discretized form of the PBE (1), suitable for formulating transformation matrices.
