*2.3. Breakage*

To ge<sup>t</sup> the equation for the pure breakage, put β(*<sup>x</sup>*, *x*) = 0 in Equation (1). Hence, the breakage process can be written as the following linear differential equation

$$\frac{\partial u(\mathbf{x},t)}{\partial t} = \int\_{\mathbf{x}}^{\infty} b(\mathbf{x},y)S(y)u(y,t)dy - S(\mathbf{x})u(\mathbf{x},t). \tag{33}$$

So, integrating Equation (33) over the cell *xi*−12 , *xi*+ 12 , can be obtained

$$\frac{d\underline{u\_i}}{dt} = B\_i - D\_{i\nu} \tag{34}$$

where

$$B\_{i} = \frac{1}{\Delta x\_{i}} \int\_{x\_{i-\frac{1}{2}}}^{x\_{i+\frac{1}{2}}} \int\_{x}^{x\_{L\_{1}+\frac{1}{2}}} b(\mathbf{x}, \boldsymbol{\varepsilon}) \mathbb{S}(\boldsymbol{\varepsilon}) u(\boldsymbol{\varepsilon}, t) d\boldsymbol{\varepsilon} d\mathbf{x} \tag{35}$$

and

$$D\_i = \frac{1}{\Delta x\_i} \int\_{\mathcal{X}\_{i-\frac{1}{2}}}^{\mathcal{X}\_{i+\frac{1}{2}}} S(\mathbf{x}) u(\mathbf{x}, t) d\mathbf{x},\tag{36}$$

where *L*1 is the number of discrete size classes.

Take the approximation of the initial data as

$$u(\mathbf{x},0) = \frac{1}{\Delta \mathbf{x}\_i} \int\_{\mathbf{x}\_{i-\frac{1}{2}}}^{\mathbf{x}\_{i+\frac{1}{2}}} u\_0(\mathbf{x})d\mathbf{x}.\tag{37}$$

*Processes* **2019**, *7*, 535

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

$$\frac{du\_i}{dt} = \frac{1}{\Delta x\_i} \sum\_{k=i}^{L\_1} \varphi\_k^b S\_k u\_k \Delta x\_k B\_{ik} - \psi\_i^d S\_i u\_i. \tag{38}$$

where the weighted parameters for birth and death ϕ*bi*and ψ*di*are given by

$$\varphi\_i^b = \frac{\mathbf{x}\_i(\mathfrak{d}(\mathbf{x}\_i) - 1)}{\sum\_{k=1}^{i-1} (\mathbf{x}\_i - \mathbf{x}\_k)B\_{ki}} \tag{39}$$

and

$$\psi\_{i}^{d} = \frac{\varphi\_{i}^{b}}{x\_{i}} \sum\_{k=1}^{i} x\_{k} B\_{ki} \quad \text{i} = 2, 3, 4, \dots, L\_{1} \; . \tag{40}$$

Take ϕ*b*1 = 0, ψ*d*1 = 0 and

$$B\_{ik} = \int\_{x\_{i-\frac{1}{2}}}^{p\_k^t} b(x, x\_k) dx \tag{41}$$

with

$$p\_k^i = \begin{cases} \mathbf{x}\_{i\nu} & k=i\\ \mathbf{x}\_{i+\frac{1}{2}\nu} & k \neq i \end{cases} \tag{42}$$

and

$$\mathcal{S}(\mathbf{x}\_{i}) = \int\_{0}^{\mathbf{x}\_{i}} b(\mathbf{x}\_{i}|\mathbf{x}\_{i})d\mathbf{x}.\tag{43}$$

Multiplying Δ*xi* on both sides of Equation (38), one obtains

$$\frac{d\Delta \mathbf{x}\_{i}\boldsymbol{u}\_{i}}{dt} = \sum\_{k=i}^{L\_{1}} \boldsymbol{\varphi}\_{k}^{b} \mathbf{S}\_{k} \boldsymbol{\mu}\_{k} \Delta \mathbf{x}\_{k} \mathbf{B}\_{ik} - \boldsymbol{\psi}\_{i}^{d} \mathbf{S}\_{i} \boldsymbol{u}\_{i} \Delta \mathbf{x}\_{i}. \tag{44}$$

Take *gi* = Δ*xiui*. Then, Equation (44) can be rewritten as

$$\frac{d\mathcal{g}\_i}{dt} = \sum\_{k=i}^{L\_1} \varphi\_k^b S\_k g\_k B\_{ik} - \psi\_i^d S\_i \mathcal{g}\_i. \tag{45}$$

To ge<sup>t</sup> the particle number at some certain time point τ, the time span [0, τ] is divided into *K* subintervals [*t*, *t* + <sup>Δ</sup>*t*], where *t*1 = 0, *tk* = τ and Δ*t* is the time step.

The particle number *g*(*<sup>t</sup>* + Δ*t*) at time point *t* + Δ*t* using transformation matrix *<sup>T</sup>*(*<sup>t</sup>*, Δ*t*) is given by

$$
\mathcal{g}(t + \Delta t) = \mathcal{g}(t) \cdot T(t, \Delta t),
\tag{46}
$$

where *g*(*t*) is a 1 × *L*1 row matrix (vectors) with elements denoted by *gi*(*t*), and *<sup>T</sup>*(*<sup>t</sup>*, Δ*t*) is a *L*1 × *L*1 matrix whose elements are denoted by *Tij*(*<sup>t</sup>*, <sup>Δ</sup>*t*). The elements of the transformation matrix can be written with the help of the discretized form of Equation (33) as follows

$$T\_{ij}(t, \Delta t) = \begin{cases} 1 + \left(\boldsymbol{\varrho}\_{\boldsymbol{i}}^{b} \boldsymbol{B}\_{\boldsymbol{i}\boldsymbol{i}} - \boldsymbol{\psi}\_{\boldsymbol{i}}^{d}\right) \mathbf{S}\_{i} \Delta t, & \boldsymbol{i} = \boldsymbol{j} \\\ \boldsymbol{S}\_{i} \boldsymbol{B}\_{\boldsymbol{j}i} \boldsymbol{\varrho}\_{\boldsymbol{i}}^{b} \Delta t, & \boldsymbol{i} > \boldsymbol{j} \\\ \boldsymbol{0}, & \boldsymbol{i} < \boldsymbol{j} \end{cases} \tag{47}$$

The scheme (Equation (46)) gives the first-order accuracy. To ge<sup>t</sup> a higher order accuracy,

$$
\lg(t + \Delta t) = \lg(t) \cdot T(t, \Delta t),
\tag{48}
$$

is proposed using the Runge–Kutta method. Here

$$\widetilde{T}(t,\Delta t) = \left[T(t,\Delta t) \cdot T\left(t + \Delta t, \frac{\Delta t}{2}\right) + \left(I - T\left(t, \frac{\Delta t}{2}\right)\right)\right].\tag{49}$$

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

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

$$
\Delta t \le \min\_i \left| \frac{2g\_i(t)}{Q(t)} \right|, \tag{50}
$$

where

$$\begin{aligned} Q(t) &= \mathcal{W}\_1(\mathcal{g}(t)) - \mathcal{W}\_2(\mathcal{g}(t)) + \mathcal{W}\_1(\mathcal{g}^\*(t)) - \mathcal{W}\_2(\mathcal{g}^\*(t)), \\ &\quad W\_1(\mathcal{g}(t)) = \sum\_{k=i}^{L\_1} \mathcal{q}\_k^b S\_k \mathcal{g}\_k(t) B\_{ik} - \psi\_i^d S\_i \mathcal{g}\_i(t), \\ &\quad \mathcal{W}\_2(\mathcal{g}(t)) = \psi\_i^d S\_i \mathcal{g}\_i(t), \\ &\quad \mathcal{g}\_i^\*(t) = \mathcal{g}\_i(t) + \Delta t \mathcal{W}\_1(\mathcal{g}(t)) - \mathcal{W}\_2(\mathcal{g}(t)). \end{aligned} \tag{51}$$
