Application of Transformation Matrices to the Solution of Population Balance Equations

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. Record Type: Published Article Submitted To: LAPSE (Living Archive for Process Systems Engineering) Citation (overall record, always the latest version): LAPSE:2019.1120 Citation (this specific file, latest version): LAPSE:2019.1120-1 Citation (this specific file, this version): LAPSE:2019.1120-1v1 DOI of Published Version: https://doi.org/10.3390/pr7080535 License: Creative Commons Attribution 4.0 International (CC BY 4.0) Powered by TCPDF (www.tcpdf.org)


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.

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.

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 1 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(x, 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-differential equation, namely a population balance equation [23,24] ∂u(x, t) 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 . Q in (t) and . Q out (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(x, 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 S(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: and y 0 xb(x, y)dx = y, ∀y > 0.
Equation (2) designates that the breakage function b(x, 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 ϑ(y) and, obviously, ϑ(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][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.

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]: Here F is a model function that depends on: • X(t), Y(t)-input and output stream variables; • c(t)-control variables that describe the operating conditions of the process and are usually controlled in certain ranges; • p(t)-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. {X}, {Y}, and {H} 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(t), the control variables c(t), and the model parameters p(t) from Equation (4) together describe the influence of the process environment on the model during the simulation. c(t) and p(t) are not shown on the scheme for simplicity, whereas all X(t) at each moment of time are denoted as {X}. Emphasizing the dynamic nature of the model, dependent variables Y(t) are divided into two assemblies: {H} is a set of variables describing the internal state, or holdup, of the model at each time point, while {Y} 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 {F H } and {F Y } 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 {F H } and {F Y } for a particular simulation.

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]: Here is a model function that depends on:  ( ), ( )-input and output stream variables;  ( )-control variables that describe the operating conditions of the process and are usually controlled in certain ranges;  ( )-model parameters, which are customizable settings of the model itself; they may or may not change over time;  -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 (4) is also divided into the two sets { } and { } 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 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 { } and { } for a particular simulation. According to Dosta [6], there are two strategies to calculate holdup and outlet streams for steadystate 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. 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: 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: 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: {F H } and {F Y } 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.
for dynamic units or 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 R N . The model will not work, since it requires more information about distributed parameters than is available in the flowsheet.
(3) R M ⊂ 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 {F H } and {F Y } 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: Processes 2019, 7, 535 6 of 29

•
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 θ of material transition between all classes: where T is the transformation matrix. Define I N as a set of indices to address any value in the N-dimensional space, so that The dimension of the transformation matrix T is twice the number of input dimensions taken into account. Each element T {I N },{J N } of T denotes a fraction of material moving from the I N -th cell of the input distribution to the J N -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 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 Nand 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 Consider the steady-state case (13) with any X ∈ X(t) and Y ∈ Y(t) . For a simple 1D case with N = 1 and R M = R N , the operation of applying the transformation matrix can be written as Processes 2019, 7, 535 7 of 29 Then, similarly for the general case with R M = R N : 15) or using Equation (11): Then for the case R M ⊂ R N , Equation (16) becomes the following form: For example, for a three-dimensional case, if N is defined for {L 1 , L 2 , L 3 } and M is defined for {L 2 }, Equation (17) can be written as 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 ⊂ 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 θ (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.

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(x) = 0 and b(x, y) = 0 in Equation (1), the general agglomeration PBE for batch process can be written as Define the initial approximation u 0 (x) as Processes 2019, 7, 535 8 of 29 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 Integrating Equation (19) over each cell , and introducing weights as given by Kumar et al. [44], the discretized form can be obtained where w b jk and w d ij are the weights responsible for preservation of number and conservation of mass, respectively. These weights are defined as and Here l ij is the symmetric index of the cell, where the agglomerate particle x i + x j falls. Rewrite Equation (22) as Take g i = ∆x i u i , then Equation (25) gives To get the particle number at some certain time point τ, the time span [0, τ] is divided into K subintervals [t, t + ∆t], where t 1 = 0, t K = τ and ∆t is the time step. Now the particle number g(t + ∆t) at time point t + ∆t using transformation matrix T(t, ∆t) is given by Here g(t) is a 1 × L 1 row matrix with i-th element g i (t), i = 1, 2, 3, . . . , L 1 at time t, and T(t, ∆t) is a L 1 × L 1 transformation matrix, calculated at time t to obtain the distribution after time step ∆t. The elements of the transformation matrix are given by Processes 2019, 7, 535 9 of 29 To improve the accuracy of the scheme (27), the Runge-Kutta method is proposed: where 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 where

Breakage
To get the equation for the pure breakage, put β(x, x ) = 0 in Equation (1). Hence, the breakage process can be written as the following linear differential equation So, integrating Equation (33) over the cell where and where L 1 is the number of discrete size classes. Take the approximation of the initial data as Proceeding in the same way as Kumar et al. [43], the discretized form with two weight functions can be obtained as where the weighted parameters for birth and death ϕ b i and ψ d i are given by and Take and Multiplying ∆x i on both sides of Equation (38), one obtains Take g i = ∆x i u i . Then, Equation (44) can be rewritten as To get the particle number at some certain time point τ, the time span [0, τ] is divided into K subintervals [t, t + ∆t], where t 1 = 0, t k = τ and ∆t is the time step.
The particle number g(t + ∆t) at time point t + ∆t using transformation matrix T(t, ∆t) is given by where g(t) is a 1 × L 1 row matrix (vectors) with elements denoted by g i (t), and T(t, ∆t) is a L 1 × L 1 matrix whose elements are denoted by T ij (t, ∆t). The elements of the transformation matrix can be written with the help of the discretized form of Equation (33) as follows The scheme (Equation (46)) gives the first-order accuracy. To get a higher order accuracy, is proposed using the Runge-Kutta method. Here 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 (0) for all i. However, the solution g i (t) may become negative. In order to ensure positivity, a condition for the time step must be defined: where

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.

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 + 3 mm 3 )/2 = 2.5 mm 3 ; 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 -dimensional distribution represented by hierarchy levels, where is the number of classes in -th dimension and and are relative mass fractions in the -th class of the -th level in the initial and the transformed distribution, respectively (Figure 3).  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 + 3 mm 3 )/2 = 2.5 mm 3 ; 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.
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 p q and Y p q are relative mass fractions in the q-th class of the p-th level in the initial and the transformed distribution, respectively (Figure 3). 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 -dimensional distribution represented by hierarchy levels, where is the number of classes in -th dimension and and are relative mass fractions in the -th class of the -th level in the initial and the transformed distribution, respectively (Figure 3). If the transformation matrix is calculated for -dimensional parameter space, its application will consist of three stages: 1. Apply to the lowest -th hierarchy level of input distribution, using Equation (17) in the form of If the transformation matrix T is calculated for M-dimensional parameter space, its application will consist of three stages: 1.
Apply T to the lowest M-th hierarchy level of input distribution, using Equation (17) in the form of . . .
For example, for M = 2, the second level of the hierarchy should be calculated as 2.
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 which corresponds to Equation (14).

3.
Apply T to calculate the remaining levels K below M as For example, if M = 1 and K = 2: 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.

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.

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 . m out (t) equals to the mass flow entering it . m in (t), and the distribution of the material at the output Y(t) equals to the distribution in the holdup H(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: Here β(x, 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 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 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 (src1, j1) and (src2, j2) 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 src1, j1 T src1,dst )· c src2, j2 ·T src2,dst whereas the class index j for equidistant grid can be obtained as j = 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 .

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: where S(x) 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: Here B(x, y) denotes the mass fraction of particles with the size of y, which after breakage get 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.

Screen
The screen apparatus is a steady-state unit, which is described in terms of the grade efficiency G(x), calculated according to the model of Plitt [61] as (63) G(x) 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; α 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 where . 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).

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

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 −3 mm 3 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. 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 −3 mm 3 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 and standard deviations shown in Table 1. 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.  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.
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 * (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.  Both particle sizes and API concentrations are given in terms of the Gaussian normal distribution with mean values µ and standard deviations σ shown in Table 1. 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.  Figures 6 and 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.
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 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. 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 β * (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. Table 2. Material transitions during agglomeration.
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.

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.

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. The mill is initially filled with 50 kg of the same blend (Figure 9). The particle size is described between 0 to 4 × 10 −3 mm 3 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. Table 3. Model parameters used for breakage process.  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. The mill is initially filled with 50 kg of the same blend ( Figure 9). 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.     The particle size is described between 0 to 4 × 10 −3 mm 3 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. Figures 10 and 11 show the distribution of the product after this time point.   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.

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 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.

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. Initial distributions of material are given by Gaussian function.  All model parameters used to simulate the process from Figure 12 can be found in Table 4. 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.  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.

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 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.

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. Q out distribution of particles by size at the input and at the output of a unit, respectively r design or structural parameters of a unit R maximum size of particles in the considered interval R N parameter space formed by N distributed parameters s 0 size-independent selection rate factor S(y) breakage rate of a particle of size y S i breakage rate of a particle from size-class i t time T(t, ∆t) transformation matrix, calculated at time t to obtain the distribution after time step ∆t T ij particular entry of a transformation matrix with specific indices T H , T Y transformation matrices to calculate holdups and output of a unit, respectively u(x) mass fraction of particles of size x u 0 (x) mass fraction of particles of size x at the initial moment of time u i mass fraction of particles from size-class i w b , w d weighting parameters for agglomeration, for birth and death of particles, respectively x, x particle sizes x a,min ,x a,max critical sizes of particles for agglomeration x b,min ,x b,max critical sizes of particles for breakage x cut cut size of the screen unit x i size of particles in a discrete class i X unit's input variables X p q relative mass fraction, stored in the q-th class of the p-th hierarchical level in the initial distribution y particle size y minimum fragment size of the Vogel breakage function Y unit's output variables Y p q relative mass fraction, stored in the q-th class of the p-th hierarchical level in the transformed distribution α separation sharpness of the screen unit β(x, x ) agglomeration rate for particles of sizes x and x β ij agglomeration rate for particles from size-classes i and j β 0 size-independent agglomeration rate constant β * size-dependent agglomeration rate constant ∆x i length of size-class i ϑ(y) total number of particles of size y appearing after breakage θ H , θ Y model functions of a unit in terms of the laws of material transition µ conc mean value of the normal distribution function describing API concentration µ size mean value of the normal distribution function describing particle sizes σ conc standard deviation of the normal distribution function describing API concentration σ size standard deviation of the normal distribution function describing particle sizes ϕ b , ψ d weighting parameters for breakage, for birth and death of particles, respectively operation of applying the transformation matrix I N set of indices to address any value in the N-dimensional parameter space