*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(t) = F(t, X(t), c(t), p(t), r). \tag{4}$$

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


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 *<sup>X</sup>*(*t*), the control variables *<sup>c</sup>*(*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 {*FH*} and {*FY*} 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 {*FH*} and {*FY*} for a particular simulation.

**Figure 1.** Simplified scheme of embedding a model into a flowsheet.

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


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{cases} \begin{aligned} H(t + \Delta t) &= F\_H(t + \Delta t, \{X(t + \Delta t)\}, \{H(t)\}, \{Y(t)\})\\ \end{aligned} \end{cases} \begin{aligned} \begin{aligned} \begin{aligned} H(t) &\downarrow (Y(t)), \{Y(t)\} \end{aligned} \end{aligned} \end{cases} \tag{5}$$

where *FH* and *FY* 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(t) = F\_Y(t, \, |X(t)\rangle). \tag{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 *Ld*. Then, if the solid material is distributed over *N* property coordinates, the total number of state variables is defined as:

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

{*FH*} and {*FY*} are usually defined to be fully flexible related to the number of discrete classes (*Ld*) 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 {*FH*}, {*FY*}, and the holdups {*H*} are defined in the *M*-dimensional parameter space. Since the output distribution {*Y*} is calculated based on the internal states {*H*} and the model functions {*FY*}, it will also be defined in the *M*–parameter space:

$$\{X\} \in \mathbb{R}^N$$

$$\{F\_H\} \lrcorner \{F\_Y\} \in \mathbb{R}^M \to \{H\} \lrcorner \{Y\} \in \mathbb{R}^M$$

for dynamic units or

$$\{X\} \in \mathbb{R}^N$$

$$\{F\_Y\} \in \mathbb{R}^M \to \{Y\} \in \mathbb{R}^M$$

for steady-state units.

> Then three cases are possible:


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


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:

$$\begin{cases} \begin{array}{l} T\_H(t) = \theta\_H[F\_H(t)] \\ T\_Y(t) = \theta\_Y[F\_Y(t)] \end{array} \end{cases} \tag{10}$$

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

$$\{I^N\} = i\_1 i\_2 \dots i\_N, \quad i\_d \in [1:L\_d], \quad d \in [1:N]. \tag{11}$$

The dimension of the transformation matrix *T* is twice the number of input dimensions taken into account. Each element *<sup>T</sup>*{*I<sup>N</sup>*},{*J<sup>N</sup>*} of *T* denotes a fraction of material moving from the *I <sup>N</sup>*–th cell of the input distribution to the *J <sup>N</sup>*–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{cases} H(t + \Delta t) = T\_H(t + \Delta t) \otimes H(t) \\ \quad \mathcal{Y}(t + \Delta t) = T\_{\mathcal{Y}}(t + \Delta t) \otimes H(t + \Delta t) \\ \quad \{T\_H\}, \{T\_{\mathcal{Y}}\} \in \mathbb{R}^M \to \{H\}, \{\mathcal{Y}\} \in \mathbb{R}^N \end{cases} \tag{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{aligned} Y(t) &= T\_Y(t) \otimes X(t) \\ \{T\_Y\} \in \mathbb{R}^M &\to \{Y\} \in \mathbb{R}^N. \end{aligned} \tag{13}$$

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

$$\forall j\_1 \in [1:L\_1]: Y\_{j\_1} = \sum\_{i\_1=1}^{L\_1} X\_{i\_1} \cdot T\_{i\_1, j\_1}.\tag{14}$$

Then, similarly for the general case with *R<sup>M</sup>* = *RN*:

$$\forall j\_1 \in [1:L\_1], \forall j\_2 \in [1:L\_2], \dots, \forall j\_N \in [1:L\_N]:$$

$$\forall j\_{j\_1 j\_2 \dots j\_N} = \sum\_{i\_1=1}^{L\_1} \sum\_{i\_2=1}^{L\_2} \dots \sum\_{i\_N=1}^{L\_N} X\_{i\_1 i\_2 \dots i\_N} \cdot T\_{i\_1 i\_2 \dots i\_N, j\_1 j\_2 \dots j\_N} \tag{15}$$

or using Equation (11):

$$\mathcal{N}\{\mathbf{J}^{\rm N}\}: \mathcal{Y}\_{\{\mathbf{J}^{\rm N}\}} = \sum\_{\{\mathbf{I}^{\rm N}\}} X\_{\{\mathbf{I}^{\rm N}\}} \Delta T\_{\{\mathbf{I}^{\rm N}\}, \{\mathbf{J}^{\rm N}\}}.\tag{16}$$

Then for the case *R<sup>M</sup>* ⊂ *<sup>R</sup>N*, Equation (16) becomes the following form:

$$\mathsf{V}\{\mathsf{J}^{\mathsf{N}}\}:\mathsf{Y}\_{\{\mathsf{J}^{\mathsf{N}}\}}=\sum\_{\{\mathsf{I}^{\mathsf{N}}\}}\mathsf{X}\_{\{\mathsf{I}^{\mathsf{N}}\}}\cdot T\_{\{\mathsf{I}^{\mathsf{M}}\},\{\mathsf{J}^{\mathsf{M}}\}}.\tag{17}$$

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

$$\forall j\_1 \in [1:L\_1], \forall j\_2 \in [1:L\_2], \forall j\_3 \in [1:L\_3]: Y\_{j\_1 j\_2 j\_3} = \sum\_{i\_1=1}^{L\_1} \sum\_{i\_2=1}^{L\_2} \sum\_{i\_3=1}^{L\_3} X\_{i\_1 i\_2 i\_3} \cdot T\_{i\_2, j\_2}.\tag{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<sup>M</sup>* ⊂ *RN*. 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.
