Application of Transformation Matrices to the Solution of Population Balance Equations

Abstract

The development of algorithms and methods for modelling flowsheets in the field of granular materials has a number of challenges. The difficulties are mainly related to the inhomogeneity of solid materials, requiring a description of granular materials using distributed parameters. To overcome some of these problems, an approach with transformation matrices can be used. This allows one to quantitatively describe the material transitions between different classes in a multidimensional distributed set of parameters, making it possible to properly handle dependent distributions. This contribution proposes a new method for formulating transformation matrices using population balance equations (PBE) for agglomeration and milling processes. The finite volume method for spatial discretization and the second-order Runge–Kutta method were used to obtain the complete discretized form of the PBE and to calculate the transformation matrices. The proposed method was implemented in the flowsheet modelling framework Dyssol to demonstrate and prove its applicability. Hence, it was revealed that this new approach allows the modelling of complex processes involving materials described by several interconnected distributed parameters, correctly taking into consideration their interdependency.

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:

• Aspen Plus [10], which is actively used in chemical industry as well as for simulation of polymers, minerals, metals etc.
• gPROMS FormulatedProducts [11], as a part of gPROMS platform, especially designed to investigate solid phase processes.
• JKSimMet [12], which is a flowsheeting software for simulation of comminution and classification circuits in the mineral processing industry.
• The HSC Sim [13] module of HSC Chemistry software intended for modelling various processes in chemistry, metallurgy, and mineralogy.
• CHEMCAD [14], which was developed to simulate chemical processes with limited consideration of the solid phase.
• Dyssol [7], a flowsheet simulation system designed to simulate complex dynamic processes in solids processing technology, developed within a research collaboration founded by the German Research Foundation.

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¹ 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 $u\left(x,t\right)\ge 0$ in these processes, at time $t\ge 0$, for particles of size $x\ge 0$ can be represented by a special type of nonlinear integro-differential equation, namely a population balance equation [23,24]

$$\begin{gathered} \frac{\partial u\left(x,t\right)}{\partial t}=\frac{1}{2}{{\displaystyle \int }}_0^x\beta \left(x^{\prime },x-x\prime \right)u\left(x^{\prime },t\right)u\left(x-x^{\prime },t\right)dx\prime -{{\displaystyle \int }}_0^{\infty }\beta \left(x,x^{\prime }\right)u\left(x,t\right)u\left(x^{\prime },t\right)d{x}^{\prime } \\ +{{\displaystyle \int }}_x^{\infty }b\left(x,y\right)S\left(y\right)u\left(y,t\right)dy-S\left(x\right)u\left(x,t\right)+{\dot{Q}}_{in}\left(t\right)-{\dot{Q}}_{out}\left(t\right). \end{gathered}$$

(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 ${\dot{Q}}_{in}\left(t\right)$ and ${\dot{Q}}_{out}\left(t\right)$ denote the input and the output in a continuous process, accordingly. The function $\beta \left(x^{\prime },x\right)$ is the agglomeration kernel, which represents the rate at which a particle of size $x$ agglomerates with a particle of size $x\prime$. The breakage function $b\left(x,y\right)$ in the third term represents the fraction of fragments of size $x$, which appear when a particle of size $y$ breaks. The selection function $S\left(x\right)$ 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{array}{ll}b\left(x,y\right)=0,\hfill & \forall x\ge y\hfill \\ {{\displaystyle \int }}_0^{y}b\left(x,y\right)dx=\vartheta \left(y\right),\hfill & \forall y>0\hfill \end{array}$$

(2)

and

$$\begin{array}{cc}\underset{0}{\overset{y}{{\displaystyle \int }}}xb\left(x,y\right)dx=y,& \forall y>0.\end{array}$$

(3)

Equation (2) designates that the breakage function $b\left(x,y\right)$ 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 $\vartheta \left(y\right)$ and, obviously, $\vartheta \left(y\right)\ge 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,26,27,28], the fixed pivot technique [29,30], the Monte Carlo simulation method [31,32,33,34], the finite element method [35,36,37], approaches based on a separable approximation of the kernel function with the fast Fourier transformation [38,39,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 get the full discretized form of the PBE (1), suitable for formulating transformation matrices.

2. Mathematical Formulations

2.1. Transformation Matrices

With regard to the flowsheet simulation, mathematical models that describe time-dependent changes occurring at individual stages in the process can be represented by the following general equation [45]:

$$Y\left(t\right)=F\left(t,X\left(t\right),c\left(t\right),p\left(t\right),r\right).$$

(4)

Here $F$ is a model function that depends on:

• $X\left(t\right)$, $Y\left(t\right)$—input and output stream variables;
• $c\left(t\right)$—control variables that describe the operating conditions of the process and are usually controlled in certain ranges;
• $p\left(t\right)$—model parameters, which are customizable settings of the model itself; they may or may not change over time;
• $r$—design variables that represent structural features of apparatuses, such as physical dimensions and shapes, and usually do not change during the simulation.

Considering this, the information flow between a flowsheet and a dynamic model can be simplistically represented by the scheme shown in Figure 1. $\left\{X\right\}$, $\left\{Y\right\}$, and $\left\{H\right\}$ are sets of time-dependent state variables, which describe inlet, outlet, and holdup of the model, respectively. They are composed of distributed and concentrated properties of bulk material, such as size, humidity, or composition of particles. The input variables $X\left(t\right)$, the control variables $c\left(t\right)$, and the model parameters $p\left(t\right)$ from Equation (4) together describe the influence of the process environment on the model during the simulation. $c\left(t\right)$ and $p\left(t\right)$ are not shown on the scheme for simplicity, whereas all $X\left(t\right)$ at each moment of time are denoted as $\left\{X\right\}$. Emphasizing the dynamic nature of the model, dependent variables $Y\left(t\right)$ are divided into two assemblies: $\left\{H\right\}$ is a set of variables describing the internal state, or holdup, of the model at each time point, while $\left\{Y\right\}$ is a set of output parameters which describe material leaving the model. The model function $F$ from Equation (4) is also divided into the two sets $\left\{F_H\right\}$ and $\left\{F_Y\right\}$ to describe changes in material parameters occurring in the holdup and in the outlet streams of the model, respectively. In contrast to all input time-dependent parameters, design variables $r$ do not change during the whole simulation. Therefore, they can be considered as model constants for each individual process and, thereby, as internal parameters of $\left\{F_H\right\}$ and $\left\{F_Y\right\}$ for a particular simulation.

According to Dosta [6], there are two strategies to calculate holdup and outlet streams for steady-state and dynamic models:

(1)

Explicitly, by direct calculation of all state variables in holdups and outlet streams.

(2)

Implicitly, through the application of movement (transformation) matrices, which describe the transfer of material between discrete classes in multidimensional parameter space.

The explicit approach is most commonly used in flowsheet software. In the case of an explicit calculation of a dynamic model using the Euler integration scheme, Equation (4) can be reformulated for models, discretized with respect to time, as:

$$\{\begin{array}{l}H\left(t+\Delta t\right)=F_H\left(t+\Delta t,\left\{X\left(t+\Delta t\right)\right\},\left\{H\left(t\right)\right\},\left\{Y\left(t\right)\right\}\right)\hfill \\ Y\left(t+\Delta t\right)=F_Y\left(t+\Delta t,\left\{H\left(t+\Delta t\right)\right\}\right)\hfill \end{array},$$

(5)

where $F_H$ and $F_Y$ are functions to describe changes occurring in holdup and outlet streams. Since steady-state models do not include holdups, their models have a simplified form:

$$Y\left(t\right)=F_Y\left(t,\left\{X\left(t\right)\right\}\right).$$

(6)

To describe a continuous distribution of material through a specific property, a discretized representation is used. Each property $d$ is described by a certain number of discrete classes $L_d$. Then, if the solid material is distributed over $N$ property coordinates, the total number of state variables is defined as:

$$L={\displaystyle \prod }_{d=1}^N{L}_d.$$

(7)

$\left\{F_H\right\}$ and $\left\{F_Y\right\}$ are usually defined to be fully flexible related to the number of discrete classes ($L_d$) for each property coordinate. Moreover, to make the simulation system generally applicable, the number of the property coordinates $N$ and their types should also be flexible. This means that the final set of parameters may not be known during model development. However, if the model is designed applying the explicit approach, the number of considered distributed properties $M$ is always fixed. That means that for the general case the model functions $\left\{F_H\right\}$, $\left\{F_Y\right\}$, and the holdups $\left\{H\right\}$ are defined in the $M$-dimensional parameter space. Since the output distribution $\left\{Y\right\}$ is calculated based on the internal states $\left\{H\right\}$ and the model functions $\left\{F_Y\right\}$, it will also be defined in the $M$–parameter space:

$$\begin{gathered} \left\{X\right\}\in R^N \\ \left\{F_H\right\},\left\{F_Y\right\}\in R^M\to \left\{H\right\},\left\{Y\right\}\in R^M \end{gathered}$$

(8)

for dynamic units or

$$\begin{gathered} \left\{X\right\}\in R^N \\ \left\{F_Y\right\}\in R^M\to \left\{Y\right\}\in R^M \end{gathered}$$

(9)

for steady-state units.

Then three cases are possible:

(1)

$R^M=R^N$. The model considers in its equations all the distributed parameters given in the inlet streams. In this case, the explicit solution works well, since all distributed parameters can be calculated directly.

(2)

$R^M\notin R^N$. The model will not work, since it requires more information about distributed parameters than is available in the flowsheet.

(3)

$R^M\subset R^N$. The explicit scheme will provide proper results in the $M$-dimensional parameter space, but will fail to deliver correct values for $R^N/R^M$, since $\left\{F_H\right\}$ and $\left\{F_Y\right\}$ cannot be properly applied for other dimensions, beyond those for which they were defined.

In the third case, all the disadvantages and limitations of the explicit approach become apparent:

• The whole set of possible distributed parameters must be known during the development of the model and they all should be considered in its equations.
• If the number or composition of the distributed parameters alters in an individual simulation, the model itself should track these changes and react to them properly.
• The model must ensure setting all distributions defined at the input to its holdup and output, even those that are not explicitly considered in the model.
• Considering only those distributed parameters that are necessary for the model leads to the loss of information about the remaining ones, despite the fact that there may be enough information to calculate them.

Thus, the use of the method based on transformation matrices is an option for simulation systems with an unlimited number of distributions, since it allows treating all variables correctly, regardless of the number and type of defined distributed parameters. In this case, the model functions are used not directly, but to derive the laws $\theta$ of material transition between all classes:

$$\{\begin{array}{l}{T}_H\left(t\right)={\theta }_H\left[F_H\left(t\right)\right]\hfill \\ T_Y\left(t\right)={\theta }_Y\left[F_Y\left(t\right)\right]\hfill \end{array},$$

(10)

where $T$ is the transformation matrix.

Define $\left\{I^N\right\}$ as a set of indices to address any value in the $N$-dimensional space, so that

$$\begin{array}{ccc}\left\{I^N\right\}=i_1{i}_2\dots i_N,& i_d\in \left[1:L_d\right],& d\in \left[1:N\right].\end{array}$$

(11)

The dimension of the transformation matrix $T$ is twice the number of input dimensions taken into account. Each element $T_{\left\{I^N\right\},\left\{J^N\right\}}$ of $T$ denotes a fraction of material moving from the $\left\{I^N\right\}$–th cell of the input distribution to the $\left\{J^N\right\}$–th cell of the output distribution. The output values are obtained by applying the transformation matrix to the corresponding input. For example, for a dynamic unit, at each time step, the holdup is calculated as the transformation of the previous state of the unit while considering the inlet and the outlet streams. The output, in turn, is calculated by applying the transformation laws to the holdup. Accordingly, Equation (5) takes the form

$$\begin{gathered} \{\begin{array}{l}H\left(t+\Delta t\right)=T_H\left(t+\Delta t\right)⭙H\left(t\right)\hfill \\ Y\left(t+\Delta t\right)=T_Y\left(t+\Delta t\right)⭙H\left(t+\Delta t\right)\hfill \end{array} \\ \left\{T_H\right\},\left\{T_Y\right\}\in R^M\to \left\{H\right\},\left\{Y\right\}\in R^N, \end{gathered}$$

(12)

where $⭙$ denotes the operation of applying the transformation laws. Due to the method of calculating and applying the transformation matrix, this method provides the correct conversion between the $N$– and $M$–parameter spaces.

Application of transformation laws for steady-state units is a direct transformation of input variables into output variables, so Equation (6) becomes

$$\begin{gathered} Y\left(t\right)=T_Y\left(t\right)⭙X\left(t\right) \\ \left\{T_Y\right\}\in R^M\to \left\{Y\right\}\in R^N. \end{gathered}$$

(13)

Consider the steady-state case (13) with any $X\in \left\{X\left(t\right)\right\}$ and $Y\in \left\{Y\left(t\right)\right\}$. For a simple 1D case with $N=1$ and $R^M=R^N$, the operation of applying the transformation matrix $⭙$ can be written as

$$\forall j_1\in \left[1:L_1\right]:Y_{j_1}={\displaystyle \sum }_{i_1=1}^{L_1}{X}_{i_1}·T_{i_1,j_1}.$$

(14)

Then, similarly for the general case with $R^M=R^N$:

$$\begin{gathered} \forall j_1\in \left[1:L_1\right],\forall j_2\in \left[1:L_2\right],\dots ,\forall j_N\in \left[1:L_N\right]: \\ {Y}_{j_1{j}_2\dots j_N}={\displaystyle \sum }_{i_1=1}^{L_1}{\displaystyle \sum }_{i_2=1}^{L_2}\dots {\displaystyle \sum }_{i_N=1}^{L_N}{X}_{i_1{i}_2\dots i_N}·T_{i_1{i}_2\dots i_N,j_1{j}_2\dots j_N} \end{gathered}$$

(15)

or using Equation (11):

$$\forall \left\{J^N\right\}:Y_{\left\{J^N\right\}}={\displaystyle \sum }_{\left\{I^N\right\}}{X}_{\left\{I^N\right\}}\Delta T_{\left\{I^N\right\},\left\{J^N\right\}}.$$

(16)

Then for the case $R^M\subset R^N$, Equation (16) becomes the following form:

$$\forall \left\{J^N\right\}:Y_{\left\{J^N\right\}}={\displaystyle \sum }_{\left\{I^N\right\}}{X}_{\left\{I^N\right\}}·T_{\left\{I^M\right\},\left\{J^M\right\}}.$$

(17)

For example, for a three-dimensional case, if $N$ is defined for $\left\{L_1,L_2,L_3 \right\}$ and $M$ is defined for $\left\{L_2\right\}$, Equation (17) can be written as

$$\forall j_1\in \left[1:L_1\right],\forall j_2\in \left[1:L_2\right],\forall j_3\in \left[1:L_3\right]:Y_{j_1{j}_2{j}_3}={\displaystyle \sum }_{i_1=1}^{L_1}{\displaystyle \sum }_{i_2=1}^{L_2}{\displaystyle \sum }_{i_3=1}^{L_3}{X}_{i_1{i}_2{i}_3}·T_{i_2,j_2}.$$

(18)

In this way, the transformation matrix calculated for the $M$-dimensional unit can be applied to calculate the $N$-dimensional output from the $N$-dimensional input, even if $R^M\subset R^N$. Equations (14)–(18) can be similarly derived for dynamic units (Equation (12)).

Usually, models in solids processing technology are described using equations that directly compute output distributions. Therefore, obtaining the transformation laws $\theta$ (Equation (10)) is an additional and nontrivial task that may require significant reformulations of existing models. This paper shows how this can be accomplished for agglomeration and breakage processes.

2.2. Agglomeration

The finite volume method developed by Kumar et al. [44] is used to get the discretized form of the pure agglomeration population balance equation. Taking $S\left(x\right)=0$ and $b\left(x,y\right)=0$ in Equation (1), the general agglomeration PBE for batch process can be written as

$$\frac{\partial u\left(x,t\right)}{\partial t}=\frac{1}{2}{{\displaystyle \int }}_0^x\beta \left(x^{\prime },x-x\prime \right)u\left(x^{\prime },t\right)u\left(x-x^{\prime },t\right)dx\prime -{{\displaystyle \int }}_0^{\infty }\beta \left(x,x^{\prime }\right)u\left(x,t\right)u\left(x^{\prime },t\right)d{x}^{\prime }$$

(19)

with initial data $u\left(x,0\right)=u_0\left(x\right)\ge 0$.

Take the whole domain as $\left[0,R\right]$. Divide the domain into $L_1$ cells with boundaries $x_{i-\frac{1}{2}}$ and $x_{i+\frac{1}{2}}$ for $i=1, 2,\dots ,L_1$ where $x_{\frac{1}{2}}=0$ and $x_{L_1+\frac{1}{2}}=R$. Take the representative of the $i$–th cell $\left[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}\right]$ as $x_i=\left(\frac{x_{i+\frac{1}{2}}+x_{i-\frac{1}{2}}}{2}\right)$ and cell length $\Delta x_i=x_{i+\frac{1}{2}}-x_{i-\frac{1}{2}}$.

Define the initial approximation $u_0\left(x\right)$ as

$$u\left(0,x\right)=\frac{1}{\Delta x_i}{{\displaystyle \int }}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_0\left(x\right)dx.$$

(20)

To get the discretized form of the agglomeration Equation (19), calculate the total birth and death in each cell. For that, collect all those particles from the lower cell which is going to a particular higher cell. This is done by defining some set of indices. Here, the index set is calculated as follows

$${\mathbb{I}}^i=\left\{\left(j,k\right)\in L_1\times L_1:x_{i-\frac{1}{2}}<\left(x_j+x_k\right)\le x_{i+\frac{1}{2}}\right\}.$$

(21)

Integrating Equation (19) over each cell $\left[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}\right]$, and introducing weights as given by Kumar et al. [44], the discretized form can be obtained

$$\frac{d{u}_i}{dt}=\frac{1}{2}{\displaystyle \sum }_{\left(j,k\right)\in {\mathbb{I}}^i}{\beta }_{jk}{u}_j{u}_k\frac{\Delta x_j\Delta x_k}{\Delta x_i}{w}_{jk}^b-{\displaystyle \sum }_{j=1}^{L_1}{\beta }_{ij}{u}_i{u}_j\Delta x_j{w}_{ij}^d.$$

(22)

where $w_{jk}^b$ and $w_{ij}^d$ are the weights responsible for preservation of number and conservation of mass, respectively. These weights are defined as

$$w_{jk}^b=\{\begin{array}{cc}\frac{x_j+x_k}{2x_{l_{jk}}-\left(x_j+x_k\right)},& x_j+x_k\le R\\ 0,& x_j+x_k>R\end{array}$$

(23)

and

$$w_{ij}^d=\{\begin{array}{cc}\frac{x_{l_{ij}}}{2x_{l_{ij}}-\left(x_i+x_j\right)},& x_i+x_j\le R\\ 0,& x_i+x_j>R\end{array}.$$

(24)

Here $l_{ij}$ is the symmetric index of the cell, where the agglomerate particle $\left(x_i+x_j\right)$ falls. Rewrite Equation (22) as

$$\frac{d\mathsf{\Delta }{x}_i{u}_i}{dt}=\frac{1}{2}{\displaystyle \sum }_{\left(j,k\right)\in {\mathbb{I}}^i}{\beta }_{jk}{u}_j{u}_k\Delta x_j\Delta x_k{w}_{jk}^b-{\displaystyle \sum }_{j=1}^{L_1}{\beta }_{ij}{u}_i{u}_j\Delta x_i\Delta x_j{w}_{ij}^d.$$

(25)

Take $g_i=\mathsf{\Delta }{x}_i{u}_i$, then Equation (25) gives

$$\frac{d{g}_i}{dt}=\frac{1}{2}{\displaystyle \sum }_{\left(j,k\right)\in {\mathbb{I}}^i}{\beta }_{jk}{g}_j{g}_k{w}_{jk}^b-{\displaystyle \sum }_{j=1}^{L_1}{\beta }_{ij}{g}_i{g}_j{w}_{ij}^d.$$

(26)

To get the particle number at some certain time point $\tau$, the time span $\left[0,\tau \right]$ is divided into $K$ subintervals $\left[t,t+\Delta t\right]$, where $t_1=0$, $t_K=\tau$ and $\Delta t$ is the time step.

Now the particle number $g\left(t+\Delta t\right)$ at time point $t+\Delta t$ using transformation matrix $T\left(t,\Delta t\right)$ is given by

$$g\left(t+\Delta t\right)=g\left(t\right)·T\left(t,\Delta t\right).$$

(27)

Here $g\left(t\right)$ is a $1\times L_1$ row matrix with $i$-th element $g_i\left(t\right)$, $i=1,2,3,\dots ,L_1$ at time $t$, and $T\left(t,\Delta t\right)$ is a $L_1\times L_1$ transformation matrix, calculated at time $t$ to obtain the distribution after time step $\Delta t$. The elements of the transformation matrix are given by

$$T_{ij}\left(t,\Delta t\right)=\{\begin{array}{cc}1+\left(\frac{1}{2}{\displaystyle \sum _{\left(i,k\right)\in {\mathbb{I}}^i}}{\beta }_{ik}{g}_k\left(t\right)w_{ik}^b-{\displaystyle \sum _{j=1}^{L_1}}{\beta }_{ij}{g}_j\left(t\right)w_{ij}^d\right)\Delta t,& i=j\\ \frac{1}{2}{\displaystyle \sum _{\left(i,k\right)\in {\mathbb{I}}^i}}{\beta }_{ik}{g}_k\left(t\right)w_{ik}^b\Delta t,& i<j\\ 0,& i>j\end{array}.$$

(28)

To improve the accuracy of the scheme (27), the Runge–Kutta method is proposed:

$$g\left(t+\Delta t\right)=g\left(t\right)·\tilde{T}\left(t,\Delta t\right)$$

(29)

where

$$\tilde{T}\left(t,\Delta t\right)=\left[T\left(t,\Delta t\right)·T\left(t+\Delta t,\frac{\Delta t}{2}\right)+\left(I-T\left(t,\frac{\Delta t}{2}\right)\right)\right].$$

(30)

Here $I$ is the identity matrix. The method given by Equations (29) and (30) gives the second order accuracy.

Since the scheme is explicit, there would be some restriction on the time step to get a non-negative solution of the scheme (Equation (29)). Therefore, the Courant-Friedrichs-Lewy (CFL) condition [46] on the time step is given by

$$\Delta t\le \underset{i}{\min }\left|\frac{2g_i\left(t\right)}{Q\left(t\right)}\right|,$$

(31)

where

$$\begin{gathered} Q\left(t\right)=W_1\left(g\left(t\right)\right)-W_2\left(g\left(t\right)\right)+W_1\left(g^{*}\left(t\right)\right)-W_2\left(g^{*}\left(t\right)\right), \\ {W}_1\left(g\left(t\right)\right)=\frac{1}{2}{\displaystyle \sum }_{\left(j,k\right)\in {\mathbb{I}}^i}{\beta }_{jk}{g}_j\left(t\right)g_k\left(t\right)w_{jk}^b, \\ {W}_2\left(g\left(t\right)\right)={\displaystyle \sum }_{j=1}^{L_1}{\beta }_{ij}{g}_i\left(t\right)g_j\left(t\right)w_{ij}^d, \\ {g}_i^{*}\left(t\right)=g_i\left(t\right)+\Delta t·W_1\left(g\left(t\right)\right)-W_2\left(g\left(t\right)\right) \end{gathered}$$

(32)

2.3. Breakage

To get the equation for the pure breakage, put $\beta \left(x,x\text{'}\right)=0$ in Equation (1). Hence, the breakage process can be written as the following linear differential equation

$$\frac{\partial u\left(x,t\right)}{\partial t}={\int }_x^{\infty }b\left(x,y\right)S\left(y\right)u\left(y,t\right)dy-S\left(x\right)u\left(x,t\right).$$

(33)

So, integrating Equation (33) over the cell $\left[x_{i-\frac{1}{2}}, x_{i+\frac{1}{2}}\right]$, can be obtained

$$\frac{d{u}_i}{dt}=B_i-D_i,$$

(34)

where

$$B_i=\frac{1}{\Delta x_i}{{\displaystyle \int }}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{{\displaystyle \int }}_x^{x_{L_1+\frac{1}{2}}}b\left(x,ϵ\right)S\left(ϵ\right)u\left(ϵ,t\right)dϵdx$$

(35)

and

$$D_i=\frac{1}{\Delta x_i}{{\displaystyle \int }}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}S\left(x\right)u\left(x,t\right)dx,$$

(36)

where $L_1$ is the number of discrete size classes.

Take the approximation of the initial data as

$$u\left(x,0\right)=\frac{1}{\Delta x_i}{{\displaystyle \int }}_{x_{i-\frac{1}{2}}}^{x_{i+\frac{1}{2}}}{u}_0\left(x\right)dx.$$

(37)

Proceeding in the same way as Kumar et al. [43], the discretized form with two weight functions can be obtained as

$$\frac{d{u}_i}{dt}=\frac{1}{\Delta x_i}{\displaystyle \sum }_{k=i}^{L_1}{\phi }_k^b{S}_k{u}_k\mathsf{\Delta }{x}_k{B}_{ik}-{\psi }_i^d{S}_i{u}_i,$$

(38)

where the weighted parameters for birth and death ${\phi }_i^b$ and ${\psi }_i^d$ are given by

$${\phi }_i^b=\frac{x_i\left(\vartheta \left(x_i\right)-1\right)}{{{\displaystyle \sum }}_{k=1}^{i-1}\left(x_i-x_k\right)B_{ki}}$$

(39)

and

$$\begin{array}{cc}{\psi }_i^d=\frac{{\phi }_i^b}{x_i}{\displaystyle \sum _{k=1}^i}{x}_k{B}_{ki},& i=2,3,4,\dots ,L_1\end{array}.$$

(40)

Take ${\phi }_1^b=0$, ${\psi }_1^d=0$ and

$$B_{ik}={{\displaystyle \int }}_{x_{i-\frac{1}{2}}}^{p_k^i}b\left(x,x_k\right)dx$$

(41)

with

$$p_k^i=\{\begin{array}{cc}{x}_i,& k=i\\ x_{i+\frac{1}{2}},& k\ne i\end{array}$$

(42)

and

$$\vartheta \left(x_i\right)={{\displaystyle \int }}_0^{x_i}b\left(x,x_i\right)dx.$$

(43)

Multiplying $\mathsf{\Delta }{x}_i$ on both sides of Equation (38), one obtains

$$\frac{d\mathsf{\Delta }{x}_i{u}_i}{dt}={\displaystyle \sum }_{k=i}^{L_1}{\phi }_k^b{S}_k{u}_k\mathsf{\Delta }{x}_k{B}_{ik}-{\psi }_i^d{S}_i{u}_i\Delta x_i.$$

(44)

Take $g_i=\mathsf{\Delta }{x}_i{u}_i$. Then, Equation (44) can be rewritten as

$$\frac{d{g}_i}{dt}={\displaystyle \sum }_{k=i}^{L_1}{\phi }_k^b{S}_k{g}_k{B}_{ik}-{\psi }_i^d{S}_i{g}_i.$$

(45)

To get the particle number at some certain time point $\tau$, the time span $\left[0,\tau \right]$ is divided into $K$ subintervals $\left[t,t+\Delta t\right]$, where $t_1=0$, $t_k=\tau$ and $\Delta t$ is the time step.

The particle number $g\left(t+\Delta t\right)$ at time point $t+\Delta t$ using transformation matrix $T\left(t,\Delta t\right)$ is given by

$$g\left(t+\Delta t\right)=g\left(t\right)·T\left(t,\Delta t\right),$$

(46)

where $g\left(t\right)$ is a $1\times L_1$ row matrix (vectors) with elements denoted by $g_i\left(t\right)$, and $T\left(t,\Delta t\right)$ is a $L_1\times L_1$ matrix whose elements are denoted by $T_{ij}\left(t,\Delta t\right)$. The elements of the transformation matrix can be written with the help of the discretized form of Equation (33) as follows

$$T_{ij}\left(t,\Delta t\right)=\{\begin{array}{cc}1+\left({\phi }_i^b{B}_{ii}-{\psi }_i^d\right)S_i\Delta t,& i=j\\ S_i{B}_{ji}{\phi }_i^b\Delta t,& i>j\\ 0,& i<j\end{array}.$$

(47)

The scheme (Equation (46)) gives the first-order accuracy. To get a higher order accuracy,

$$g\left(t+\Delta t\right)=g\left(t\right)·\tilde{T}\left(t,\Delta t\right),$$

(48)

is proposed using the Runge–Kutta method. Here

$$\tilde{T}\left(t,\Delta t\right)=\left[T\left(t,\Delta t\right)·T\left(t+\Delta t,\frac{\Delta t}{2}\right)+\left(I-T\left(t,\frac{\Delta t}{2}\right)\right)\right].$$

(49)

$I$ is the identity matrix. The method given by Equations (48) and (49) gives the second order accuracy.

Consider a positive initial data $g_i\left(0\right)$ for all $i$. However, the solution $g_i\left(t\right)$ may become negative. In order to ensure positivity, a condition for the time step must be defined:

$$\Delta t\le \underset{i}{\min }\left|\frac{2g_i\left(t\right)}{Q\left(t\right)}\right|,$$

(50)

where

$$\begin{gathered} Q\left(t\right)=W_1\left(g\left(t\right)\right)-W_2\left(g\left(t\right)\right)+W_1\left(g^{*}\left(t\right)\right)-W_2\left(g^{*}\left(t\right)\right), \\ {W}_1\left(g\left(t\right)\right)={\displaystyle \sum }_{k=i}^{L_1}{\phi }_k^b{S}_k{g}_k\left(t\right)B_{ik}-{\psi }_i^d{S}_i{g}_i\left(t\right), \\ {W}_2\left(g\left(t\right)\right)={\psi }_i^d{S}_i{g}_i\left(t\right), \\ {g}_i^{*}\left(t\right)=g_i\left(t\right)+\Delta t·W_1\left(g\left(t\right)\right)-W_2\left(g\left(t\right)\right). \end{gathered}$$

(51)

3. Implementation Details

3.1. Dynamic Flowsheet Simulation in Dyssol

In this work, the Dyssol modelling framework [7] was used to perform simulations using the newly proposed methods. The system is being developed within the priority program of the German Research Foundation (DFG) SPP 1679 “Dynamic simulation of interconnected solids processes DYNSIM-FP” [47].

Dyssol is written entirely in the C++ programming language [48] using a modular architecture and an object-oriented approach [49]. This made it possible to split the entire program into separate, relatively independent modules, to facilitate their development. They are integrated into a single system using specially designed software interfaces. In particular, each model (e.g. agglomerator, mill, or screen) is an independent software library that can be created separately, and then connected to the simulation program. This greatly increases the flexibility and the extensibility of the entire system and simplifies its development and maintenance.

Due to the mathematical complexity and diversity of models in the area of solids processing technology, the calculations in Dyssol are carried out according to the sequential-modular approach [50]. Each model is calculated independently of the others, using its own most appropriate equations solvers and computational methods. The task of the simulation system itself in this case is to manage the information and material flows between the models in accordance with the flowsheet structure and to ensure convergence. The sequential-modular approach simplifies the development of new models, since they are less limited to the implementation and functionality of the modelling framework. Moreover, this method can be effectively applied even at the scale of individual process units [51].

Application of the modular approach imposes some restrictions on the flowsheet structures. In particular, schemes containing recycle flows cannot be calculated directly and require additional treatment. Before starting the simulation, all such flows must be identified and initialized with some predefined values in order to break the cycles [52]. After that, all units inside the recycle loop are calculated iteratively until convergence. To speed up this process, the initial parameters of recycles are calculated anew at each iteration using convergence methods [7,53].

Reaching convergence over a relatively large time interval can be very difficult and can become a resource-intensive task. To overcome this problem, dynamic modelling of flowsheets with recycle flows is performed using an approach based on the waveform relaxation method [54]. Here, iterative calculations are performed only for a certain short time interval [6,15], where convergence can be reached much faster. When the calculation of the current time interval is completed, the system uses data interpolation methods [7] to initialize the next time window and proceeds to its calculation.

3.2. Application of Transformation Matrices

From a mathematical point of view, the application of a full transformation matrix (Equation (16)) is a relatively simple operation and can be described as the regular multiplication of matrices. However, with a large number of matrix dimensions and a large number of classes, such a naive calculation will be very demanding in terms of computational capabilities and memory consumption. Moreover, the need to work with time-dependent parameters imposes additional restrictions on data structures. Therefore, all distributed parameters in Dyssol are stored in sparse format, represented by specially developed tree data structures [7], as is schematically shown in Figure 2 for a single time point.

Here, each hierarchy level describes the distribution of material over a specific parameter, and each entry is the mass fraction of material that has a specific combination of parameters. Since distributions often contain a lot of zero values, such a structure can drastically reduce the amount of data, which needs to be stored and processed. In addition, this data structure allows the use of relative mass fractions for each level of the hierarchy, instead of operating with absolute mass fractions. Thus, from Figure 2 it follows that there are 30% of particles with a size of (2 mm³ + 3 mm³)/2 = 2.5 mm³; of these, 50% have a sphericity of (0.5 + 0.75)/2 = 0.625; and from them, 20% have a porosity of (0.33 + 0.66)/2 = 0.5. The total percentage of particles with this combination of parameters is 0.3 × 0.5 × 0.2 = 3%.

The use of absolute mass fractions and tree structures allows for working with individual distributions without changing information in the others. So, transformation matrices can be effectively applied for each level separately and this will not affect all distributions located at higher hierarchical levels.

Take the $N$–dimensional distribution represented by $N$ hierarchy levels, where $L_d$ is the number of classes in $d$–th dimension and $X_q^p$ and $Y_q^p$ are relative mass fractions in the $q$–th class of the $p$–th level in the initial and the transformed distribution, respectively (Figure 3).

If the transformation matrix $T$ is calculated for $M$–dimensional parameter space, its application will consist of three stages:

• Apply $T$ to the lowest $M$–th hierarchy level of input distribution, using Equation (17) in the form of

  $$Y_{j_1{j}_2\dots j_{M-1}{j}_M}^M={\displaystyle \sum }_{i_1=1}^{L_1}{\displaystyle \sum }_{i_2=1}^{L_2}\dots {\displaystyle \sum }_{i_{M-1}=1}^{L_{M-1}}{\displaystyle \sum }_{i_M=1}^{L_M}{X}_{i_1{i}_2\dots i_M}^M·X_{i_1{i}_2\dots i_{M-1}}^{M-1}·\dots ·X_{i_1{i}_2}^2·X_{i_1}^1·T_{i_1{i}_2\dots i_M,j_1{j}_2\dots j_M}.$$

  (52)
  For example, for $M=2$, the second level of the hierarchy should be calculated as

  $$Y_{j_1{j}_2}^2={\displaystyle \sum }_{i_1=1}^{L_1}{\displaystyle \sum }_{i_2=1}^{L_1}{X}_{i_1{i}_2}^2·X_{i_1}^1·T_{i_1{i}_2,j_1{j}_2}$$

  (53)
• Extract the transformation laws for the previous level ($M-1$) from $T$ and apply them using Equation (52). Repeat it up to the hierarchy level 1. For example, for $M=2$, the first level will be calculated as

  $$Y_{j_1}^1={\displaystyle \sum }_{i_1=1}^{L_1}{X}_{i_1}^1·T_{i_1,j_1},$$

  (54)

  which corresponds to Equation (14).
• Apply $T$ to calculate the remaining levels $K$ below $M$ as

  $$\begin{gathered} Y_{j_1{j}_2\dots j_{K-1}{j}_K}^K=\frac{{{\displaystyle \sum }}_{i_1=1}^{L_1}{X}_{i_1{j}_2\dots j_K}^K·X_{i_1{j}_2\dots j_{K-1}}^{K-1}·\dots ·X_{i_1{j}_2}^2·X_{i_1}^1·T_{i_1,j_1}}{Y_{j_1{j}_2\dots j_{K-1}}^{K-1}·Y_{j_1{j}_2\dots j_{K-2}}^{K-2}·\dots ·Y_{j_1{j}_2}^2·Y_{j_1}^1}, \\ M<K\le N \end{gathered}$$

  (55)
  For example, if $M=1$ and $K=2$:

  $$Y_{j_1{j}_2}^2=\frac{{{\displaystyle \sum }}_{i_1=1}^{L_1}{X}_{i_1{j}_2}^2·X_{i_1}^1·T_{i_1,j_1}}{Y_{j_1}^1}.$$

  (56)

Thus, using Equations (52) and (55), as well as sorting the levels in the hierarchical data structure, one can apply the transformation matrix of any complexity to any input distribution.

Depending on the unit type, the conventional algorithm for applying transformation matrices varies:

• For steady state units: copy the distribution from the input to the output, then apply the transformation matrix to the output.
• For dynamic units: use the input distribution to calculate the holdup, apply the transformation matrix to the holdup, and then calculate the output.

In the latter case, because of time discretization, the mixing of the input material with the holdup and the application of the transformation matrix are performed sequentially.

3.3. Implementation of Units

To verify and investigate the proposed new methods, the Dyssol simulation system has been extended with several models that have been implemented using the provided program interfaces.

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 ${\dot{m}}_{out}\left(t\right)$ equals to the mass flow entering it ${\dot{m}}_{in}\left(t\right)$, and the distribution of the material at the output $Y\left(t\right)$ equals to the distribution in the holdup $H\left(t\right)$.

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{gathered} \beta \left(x,x^{\prime }\right)={\beta }_0{\beta }^{*}\left(x^{\frac{1}{3}}+{x^{\prime }}^{\frac{1}{3}}\right)\left(x^{-\frac{1}{3}}+{x^{\prime }}^{-\frac{1}{3}}\right), \\ {\beta }^{*}=\{\begin{array}{ll}1,\hfill & x+x^{\prime }\le x_{a,min} \hfill \\ 1-\frac{\left(x+x^{\prime }\right)-x_{a,min}}{x_{a,max}-x_{a,min}},\hfill & x_{a,min}<x+x^{\prime }<x_{a,max}\hfill \\ 0,\hfill & x+x^{\prime }\ge x_{a,max}\hfill \end{array}. \end{gathered}$$

(57)

Here $\beta \left(x,x^{\prime }\right)$ denotes the agglomeration rate between particles $x$ and $x^{\prime }$; ${\beta }_0$ is the size-independent agglomeration rate constant, dependent on operating conditions; ${\beta }^{*}$ is the size- dependent agglomeration rate constant; and $x_{a,min}$ and $x_{a,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 $x_{src1}$ and $x_{src2}$ aggregate into a particle of size $x_{dst}$. Then $T_{src1,dst}$ and $T_{src2,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={\displaystyle \sum }_{j=1}^{L_2}{c}_{i,j},$$

(58)

where $c_{i,j}$ is a mass fraction of granular material with size $x_i$ and value of the secondary dimension falling into class $j$. If particles from two classes $\left(src1,j1\right)$ and $\left(src2,j2\right)$ aggregate, the result lands into some class $\left(dst,j\right)$. 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_{dst,j}={\displaystyle \sum }_{j1=1}^{L_2}{\displaystyle \sum }_{j2=1}^{L_2}\frac{(c_{src1,j1}{T}_{src1,dst})·\left(c_{src2,j2}·T_{src2,dst}\right)}{{{\displaystyle \sum }}_{k=1}^{S_2}{c}_{src2,k}·T_{src2,dst}}$$

(59)

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

$$j=\frac{j1·c_{src1,j1}·T_{src1,dst}+j2·c_{src2,j2}·T_{src2,dst}}{c_{src1,j1}·T_{src1,dst}+c_{src2,j2}·T_{src2,dst}}.$$

(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⁻¹⁵.

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\left(x\right)=s_0·\{\begin{array}{ll}0, \hfill & x\le x_{b,min}\hfill \\ 1-{\left(\frac{x_{b,max }-x}{x_{b,max }-x_{b,min}}\right)}^n,\hfill & x_{b,min}<x<x_{b,max}\hfill \\ 1,\hfill & x\ge x_{b,max}\hfill \end{array},$$

(61)

where $S\left(x\right)$ is the mass fraction of particles of size $x$ that are crushed; $s_0$ introduces a size-independent selection rate factor; $x_{b,min}$ and $x_{b,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\left(x,y\right)=\{\begin{array}{ll}0.5\frac{q}{y}{\left(\frac{x}{y}\right)}^{q-2}\left(1+tanh\left(\frac{y-y\text{'}}{y\text{'}}\right)\right),\hfill & y\ge x\hfill \\ 0,\hfill & y<x\hfill \end{array}.$$

(62)

Here $B\left(x,y\right)$ denotes the mass fraction of particles with the size of $y$, which after breakage get the size less or equal to $x$; $y\text{'}$ is the minimum fragment size that can be achieved after milling; and $q$ is a power law exponent of the Vogel breakage function.

3.3.3. Screen

The screen apparatus is a steady-state unit, which is described in terms of the grade efficiency $G\left(x\right)$, calculated according to the model of Plitt [61] as

$$G\left(x\right)=1-exp\left(-0.693{\left(\frac{x}{x_{cut}}\right)}^{\alpha }\right).$$

(63)

$G\left(x\right)$ determines the mass fraction of material of size $x$, which leaves the screen through an outlet for coarse particles; $x_{cut}$ is the cut size of the classification model; $\alpha$ is the separation sharpness.

For proper consideration of the distributed parameters, each deck of the screen unit formulates two diagonal transformation matrices for the particle size distribution (PSD): for coarse ($T^c$) and for fines ($T^f$) output, so that

$$\begin{gathered} T_{i,i}^c=G\left(x_i\right) \\ {T}_{i,i}^f=1-G\left(x_i\right). \end{gathered}$$

(64)

The total mass flows of coarse ${\dot{m}}_{out}^c$ and fines ${\dot{m}}_{out}^f$ outlets can be calculated as

$$\begin{gathered} {\dot{m}}_{out}^c={\dot{m}}_{in}{\displaystyle \sum }_{i=1}^{L_1}G\left(x_i\right)·c_i \\ {\dot{m}}_{out}^f={\dot{m}}_{in}\left(1-{\displaystyle \sum }_{i=1}^{L_1}G\left(x_i\right)·c_i\right), \end{gathered}$$

(65)

where ${\dot{m}}_{in}$ is the mass flow at the inlet, $c_i$ is the mass fraction of particles of size $x_i$, and $L_1$ is the number of size classes.

Application of the transformation matrices was performed according to Equations (52) and (55).

4. Simulation Examples

All simulations were performed in the dynamic flowsheet simulation system Dyssol using the developed models, described in Section 3.3.

4.1. Agglomeration

A simple pharmaceutical process of agglomerating blends with two different concentrations of the active pharmaceutical ingredient (API) was simulated. The process structure is illustrated in Figure 4. The solid material is described using two interdependent distributed parameters:

• Particle size representing the volume of particles ranging from 0 to 4 × 10⁻³ mm³ and distributed over 100 equidistant classes; and
• the API concentration, which describes the mass content of an active ingredient from 0 to 10% divided into 500 equidistant classes.

The initial blend entering the agglomerator is a mixture of materials from two feeders (Figure 5):

• Feeder 1 supplies smaller particles with a lower API concentration;
• Feeder 2 supplies larger particles with a higher API concentration.

Both particle sizes and API concentrations are given in terms of the Gaussian normal distribution with mean values $\mu$ and standard deviations $\sigma$ shown in

Table 1. Model parameters used for the agglomeration process.

Feeder 1
Mass Flow	$\dot{m}$	0.25 kg/s
Mean value of the particle size distribution	${\mu }_{size}$	0.78 × 10−3 mm3
Standard deviation of the particle size distribution	${\sigma }_{size}$	0.8 × 10−4 mm3
Mean value of the API concentration distribution	${\mu }_{conc}$	0.025
Standard deviation of the API concentration distribution	${\sigma }_{conc}$	0.004
Feeder 2
Mass Flow	$\dot{m}$	0.25 kg/s
Mean value of the particle size distribution	${\mu }_{size}$	1.18 × 10−3 mm3
Standard deviation of the particle size distribution	${\sigma }_{size}$	0.8 × 10−4 mm3
Mean value of the API concentration distribution	${\mu }_{conc}$	0.075
Standard deviation of the API concentration distribution	${\sigma }_{conc}$	0.004
Agglomerator
Holdup mass	$m$	200 kg
Size-independent rate constant	${\beta }_0$	2 × 10−16
Minimum agglomeration size	$x_{a,min}$	0.4 × 10-3 mm3
Maximum agglomeration size	$x_{a,max}$	3.75 × 10−3 mm3

. Feeders continuously supply the material at a rate of 0.25 kg/s each. The final distribution after their mixing is shown in Figure 5. The agglomerator is initially filled with 200 kg of the same blend. All model parameters used to simulate the agglomeration process are listed in Table 1.

Figure 6 and Figure 7 show the distribution of the product stream in steady-state, which was reached after 25 min of process time. Clearly visible individual peaks, which were formed as a result of the agglomeration of particles with different sizes and different API concentrations. Since the initial blend is supplied continuously at a constant rate during the entire agglomeration process, the initial distributions A and B (Figure 7) are still distinguishable in the final mixture. Agglomeration of material B with itself gives F, which has a doubled mean size and the same API concentration as B, as expected. In turn, agglomeration of B with F (having the same concentrations) gives a peak at M, which does not alter API content and has a particle size equal to (B + F). The same is true for the material from feeder 1, agglomerating with itself according to the schema: (A + A) = C, (A + C) = E, and (A + E) + (C + C) = H.

More complex behaviour is observed if particles with different API concentrations agglomerate. In this case, the particle sizes are added up and the concentration is averaged according to Equation (59). For example, agglomeration of particles from A and B gives D. Further, G is formed as a mixture of (A + D) and (B + C). The result of the agglomeration of B and C is shifted below the concentration value of 0.05 due to the mutual influence of two factors. First, the introduction of the size-dependent growth rate ${\beta }^{*}$ (Equation (57)) leads to faster growth of smaller particles. Second, averaging by the second distribution takes into account the actual amount of material being agglomerated in selected classes (see Equations (59) and (60)).

The complete scheme of material transitions during the agglomeration, shown in Figure 7, is represented in

Table 2. Material transitions during agglomeration.

+	A	B	C	D	E	F	G
A	C	D	E	G	H	I	K
B		F	G	J	L	M	-
C			H	K	-	-	-

.

Thus, information on the change in particle size during the agglomeration process, obtained from the generated transformation matrix, made it possible to take into account and correctly process the dependent distributed parameters of the agglomerated material.

4.2. Breakage

To verify the proposed new method for solving the breakage PBE, a simplified process of milling a two-component blend was simulated. Figure 8 shows its flowsheet structure. Two feeders supply material with different particle sizes and API concentrations into the system. Their mixture, shown in Figure 9, enters the mill unit to be mixed with the holdup material and crushed afterwards.

The initial blend entering mill is a mixture of smaller particles with a higher API concentration and larger particles with a lower API concentration. Both distributed parameters are given as Gaussian normal distributions with parameters from

Table 3. Model parameters used for breakage process.

Feeder 1
Mass Flow	$\dot{m}$	0.25 kg/s
Mean value of the particle size distribution	${\mu }_{size}$	3.18 × 10−3 mm3
Standard deviation of the particle size distribution	${\sigma }_{size}$	0.8 × 10−4 mm3
Mean value of the API concentration distribution	${\mu }_{conc}$	0.025
Standard deviation of the API concentration distribution	${\sigma }_{conc}$	0.004
Feeder 2
Mass Flow	$\dot{m}$	0.25 kg/s
Mean value of the particle size distribution	${\mu }_{size}$	2.78 × 10−3 mm3
Standard deviation of the particle size distribution	${\sigma }_{size}$	0.8 × 10−4 mm3
Mean value of the API concentration distribution	${\mu }_{conc}$	0.075
Standard deviation of the API concentration distribution	${\sigma }_{conc}$	0.004
Mill
Holdup mass	$m$	50 kg
Selection rate factor	$s_0$	3 × 10−2
Minimum breakage size	$x_{b,min}$	1.5 × 10−3 mm3
Maximum breakage size	$x_{b,max}$	4 × 10−3 mm3
Selection power law exponent	$n$	3.1
Minimum fragment size	$y\prime$	0.5 × 10−3 mm3
Breakage power law exponent	$q$	5.5

. The mill is initially filled with 50 kg of the same blend (Figure 9).

The particle size is described between 0 to 4 × 10⁻³ mm³ with 100 equidistant classes, whereas the mass content of an active ingredient is distributed over 500 equidistant classes ranging from 0 to 10%.

All model parameters used to simulate breakage process are listed in Table 3.

The steady-state is reached in 25 min of the process time. Figure 10 and Figure 11 show the distribution of the product after this time point.

Due to the constant supply of material during the whole process time, initial distributions A and B (Figure 11) are still clearly visible in the product. At the same time, it is evident from Figure 10 that the reduction in particle size caused by the breakage process does not lead to changes in the content of the active component: the distribution of the API concentration in the outlet stream remains the same as in the initial distribution. Both A and B do not mix, but only change particle size gradually grinding down from A to C and from B to D.

Thus, with the help of transformation matrices, it was possible to avoid mixing the secondary parameters, despite the fact that the model of the mill unit was developed considering only the particle size distribution.

4.3. Coupled Agglomeration and Breakage

Having both models of agglomerator and mill, it is possible to simulate part of a complex pharmaceutical process with a closed circuit and an external classification of the material, as it is shown in Figure 12. New material enters the process through two feeders, supplying particles with different sizes and concentrations of API, which mixture is shown in Figure 13. This blend is mixed with the holdup inside the agglomerator, where the growth of particles occurs. The two-compartment screen unit separates material leaving the agglomerator into three fractions. Oversized particles are crushed using the mill unit, and afterwards combined with the undersized particles and sent back to the agglomerator by mixing with external material from feeders. The middle-sized agglomerates are considered as product. Initially, both agglomerator and mill are filled with the same material as it comes from feeder 1.

All model parameters used to simulate the process from Figure 12 can be found in

Table 4. Model parameters used for coupled agglomeration and breakage process.

Feeder 1
Mass Flow	$\dot{m}$	0.15 kg/s
Mean value of the particle size distribution	${\mu }_{size}$	0.78 × 10−3 mm3
Standard deviation of the particle size distribution	${\sigma }_{size}$	2.0 × 10−4 mm3
Mean value of the API concentration distribution	${\mu }_{conc}$	0.03
Standard deviation of the API concentration distribution	${\sigma }_{conc}$	0.007
Feeder 2
Mass flow	$\dot{m}$	0.35 kg/s
Mean value of the particle size distribution	${\mu }_{size}$	1.58 × 10−3 mm3
Standard deviation of the particle size distribution	${\sigma }_{size}$	3.2 × 10−4 mm3
Mean value of the API concentration distribution	${\mu }_{conc}$	0.08
Standard deviation of the API concentration distribution	${\sigma }_{conc}$	0.005
Agglomerator
Holdup mass	$m$	200 kg
Size-independent rate constant	${\beta }_0$	2 × 10-16
Minimum agglomeration size	$x_{a,min}$	0.4 × 10−3 mm3
Maximum agglomeration size	$x_{a,max}$	3.75 × 10−3 mm3
Mill
Holdup mass	$m$	50 kg
Selection rate factor	$s_0$	3 × 10-2
Minimum breakage size	$x_{b,min}$	1.5 × 10−3 mm3
Maximum breakage size	$x_{b,max}$	3.2 × 10−3 mm3
Selection power law exponent	$n$	3.1
Minimum fragment size	$y\prime$	10−3 mm3
Breakage power law exponent	$q$	5.5
Screen deck 1
Cut size	$x_{cut}$	3.2 × 10−3 mm3
Separation sharpness	$\alpha$	13
Screen deck 2
Cut size	$x_{cut}$	2.3 × 10−3 mm3
Separation sharpness	$\alpha$	13

. Initial distributions of material are given by Gaussian function.

The flowsheet was simulated until a steady state was reached, which occurred after 40 min of the process time. The distribution of the product is shown in Figure 14. There are two peaks for particles with a higher API concentration. One is mainly due to the presence of the source blend from feeder 2 (its coarse part obtained after the sieving), mixed with new particles formed in agglomerator and mill. The second peak is the result of agglomeration of particles with the same high API concentration with each other. For a low API content, only one peak is observed, since all smaller fractions (from agglomeration of the primary material from feeder 1 and milling results) were cut off by the screen unit. At the same time, a mixture of two initial blends, formed as a result of the agglomeration of materials with different API content, is present in a significant amount, but does not dominate the material with initial concentrations.

In general, using the description of granular materials with multidimensional distributions and units developed on the basis of transformation matrices, it is possible to perform simulations of complex processes involving materials distributed over a large number of dimensions. Thus, in this example, it is possible to assess the homogeneity of the material that appears at the output of the simulated process, since this parameter is crucial in pharmaceutical production.

5. Conclusions

The use of transformation matrices is a prerequisite for working with multidimensional distributed parameters of solids. This approach allows for proper handling of all parameters of granular materials, even those which are not directly handled by the model. Consequently, the number of distributed parameters in the model and in its environment may not match. At the same time, the development of each individual model becomes more complicated due to the fact that in contrast to the traditional approach, when the output distribution is explicitly calculated and set to the outlet, it is required to derive the laws of material transformation between inputs and outputs for each time step.

In this article, the applicability of transformation matrices to the description of processes of agglomeration and crushing, which are usually formulated by the population balance equations, was shown. The proposed new methods allow for extracting information about the transformation laws from PBEs and calculating the transformation matrices on that basis. The approach uses the finite volume method for spatial discretization and the second-order Runge–Kutta method for obtaining the complete discretized form of the population balance equations.

To perform the test studies, the derived methods have been implemented in the dynamic flowsheet simulation environment Dyssol. Some exemplary production processes and their parts were simulated using the derived models to prove the applicability and advantages of presented approaches for modelling of granular materials described by multidimensional distributed parameters.

As case studies have shown, the use of transformation matrices enables the correct calculation of the dependent secondary parameters, avoiding their mixing during the simulation. The proposed method can be extended with additional laws for dependent parameters, for example, for mixing of materials. This will significantly expand its application scope and will allow simulating more sophisticated processes. In general, the proposed approach allows usage of a more complex description of materials with the help of several distributed parameters, without the need to adjust each specific model, and transferring the task of their correct calculation to the modelling system. In this way, the new method of reformulating population balance equations greatly increases the possibilities for creating generally applicable systems for the simulation of solid phase processes.