1 Introduction

In the past few decades, population balance equations have been investigated and employed for modelling real life applications including depolymerization [1, 27], aerosol [3], granulation [6, 24], crystallization [15, 26] and liquid-liquid dispersions [19] related to particle technology. Mathematically, these equations are represented by an integro-partial differential equation that arises due to aggregation, fragmentation, growth and nucleation mechanisms. Among all these processes, the fragmentation event is of our interest in the current work. Fragmentation is a size reduction process during which particles of smaller size (or volume) are formed after breaking the bigger particles. Fragmentation process produces particles of varying sizes within the system, which can be tracked by the temporal change in the number density function via the population balance equation.

1.1 Fragmentation equation

The temporal change of particle number density, \(n(t, v)\ge 0\) of volume \(v \ge 0\) at time \(0\le t \le T<\infty \) in a spatially homogeneous physical system undergoing fragmentation is written as

$$\begin{aligned} \frac{\partial n(t, v)}{\partial t} = \underbrace{\int _v^\infty n(t, y)\Omega _t\left( y\right) \beta (v, y){\textrm{d}}y}_{\mathscr {B}(v)} - \underbrace{n(t, v)\Omega _t(v)}_{\mathscr {D}(v)}, \end{aligned}$$
(1.1)

and is supplemented by the initial condition (IC)

$$\begin{aligned} n(0, v) = n_0(v) \ge 0. \end{aligned}$$
(1.2)

In Eq. (1.1), \(\beta (v, y)\) is the breakage function and satisfies the following properties;

$$\begin{aligned} \int _0^y v\beta \left( v, y \right) {\textrm{d}}v = y, \quad \int _0^y \beta \left( v, y \right) {\textrm{d}}v = \zeta (y) \ge 2 \quad \text{ and }\quad \beta (v, y) = \beta (y-v, y). \end{aligned}$$
(1.3)

The first relation in the above expression (1.3) describes the mass conservation property, that is, when the particle of size y splits into smaller fragments then the total mass of those fragments must be equal to y. The second relation defines the total number of fragments \(\zeta (y)\) produced during the breakup of particles size y. The last expression, indicates the symmetric nature of the breakage function.

Additionally, \(\displaystyle \Omega _p(v,y)\) is the partial breakage kernel with an incorporated binary size distribution, which is mirror symmetrical about \(\displaystyle v = \frac{y}{2}\) and defined by the product of the total breakup rate and the daughter size distribution

$$\begin{aligned} \Omega _p(v, y) = \Omega _t\left( y\right) \beta (v, y) \quad \text{ and } \quad \Omega _t(y) = \frac{1}{2}\int _0^y \Omega _p(v, y) {\textrm{d}}v. \end{aligned}$$
(1.4)

The first term on the right-hand side of Eq. (1.1) denotes the accumulation (or birth) of particle size v in the system, whereas the second term defines the loss (or death) of particle size v from the system. For fragmentation equations, number density function is not the only important parameter to be approximated. Several integral properties of the number density function, known as moments, indicate significant real life properties which also need to be predicted with good accuracy. Therefore, the moments are defined as

$$\begin{aligned} M_i(t)=\int _{0}^{\infty }v^i n(t, v) {\textrm{d}}v, \end{aligned}$$
(1.5)

where \(M_0(t)\) (zeroth moment) and \(M_1(t)\) (first moment) provide the total number and total mass of particles in the system during a multiple fragmentation event, respectively. The focus of this study is to deal with the binary fragmentation equation which estimates the number density, zeroth and first moments with improved accuracy. In the following discussion, the zeroth and first order moments for the binary fragmentation Eq. (1.1) are calculated.

1.2 Zeroth and first order moments for the case of binary fragmentation

Consider that y breaks in two components v and \(y-v=x\) (say). Then

  1. (i)

    \(\displaystyle v\le \frac{y}{2}\), implies \(\displaystyle x \ge \frac{y}{2}\), therefore \(\displaystyle \theta \left( v\le \frac{y}{2}\right) \) is true \((=1)\), and \(\displaystyle \theta \left( x\le \frac{y}{2}\right) \) is false \((=0)\), and

  2. (ii)

    \(\displaystyle v\ge \frac{y}{2}\), implies \(\displaystyle x \le \frac{y}{2}\), therefore \(\displaystyle \theta \left( v\le \frac{y}{2}\right) \) is false \((=0)\), and \(\displaystyle \theta \left( x\le \frac{y}{2}\right) \) is true \((=1)\),

where \(\theta \) is the step function \(\displaystyle \theta (v) := \left\{ \begin{array}{ll} 0, &{}\text{ when } v \text{ is } \text{ false, } \\ 1, &{}\text{ when } v \text{ is } \text{ true. } \end{array} \right. \)

Thus the first term of the main Eq. (1.1) is redefined as

$$\begin{aligned} \mathscr {B}(v)&= \underbrace{\int _v^\infty n(t, y)\Omega _t\left( y\right) \beta (v, y)\theta \left( v\le \frac{y}{2}\right) {\textrm{d}}y}_{\mathscr {B}_1(v)} \nonumber \\&\quad + \underbrace{\int _v^\infty n(y)\int _0^{y} \Omega _t\left( y\right) \beta \left( x, y\right) \delta \left[ x- \left( y-v\right) \right] \theta \left( x\le \frac{y}{2}\right) {\textrm{d}}x{\textrm{d}}y}_{\mathscr {B}_2(v)}. \end{aligned}$$
(1.6)

To summarize, the first component \(\mathscr {B}_1(v)\) represents the accumulation of particle sizes smaller than \(\displaystyle \frac{y}{2}\), and the second component \(\mathscr {B}_2(v)\) represents the accumulation of particle sizes greater than \(\displaystyle \frac{y}{2}\), through the complementary particles x.

Mass conservation: We first focus on the term \(\mathscr {B}_1(v)\). Changing the order of integration

$$\begin{aligned} \int _0^\infty v\mathscr {B}_1(v){\textrm{d}}v&= \int _0^\infty v \int _v^\infty n(t, y)\Omega _t(y) \beta (v,y) \theta \left( v\le \frac{y}{2} \right) {\textrm{d}}y {\textrm{d}}v\\&= \int _0^\infty n(t, y)\Omega _t(y) \left[ \int _0^y v\beta (v,y) \theta \left( v\le \frac{y}{2} \right) {\textrm{d}}v\right] {\textrm{d}}y. \end{aligned}$$

Now \(\displaystyle \theta \left( v\le \frac{y}{2} \right) \) represents the accumulation of particle sizes smaller than \(\displaystyle \frac{y}{2}\), so we can write

$$\begin{aligned} \int _0^\infty v\mathscr {B}_1(v){\textrm{d}}v = \int _0^\infty n(t, y)\Omega _t(y) \left[ \int _0^{\displaystyle \frac{y}{2}} v\beta (v,y) {\textrm{d}}v\right] {\textrm{d}}y. \end{aligned}$$

We next consider \(\mathscr {B}_2(v)\), and note that it represents the accumulation of particles v which are larger than \(\displaystyle \displaystyle \frac{y}{2}\). Therefore, \(\displaystyle x=y-v\le \displaystyle \frac{y}{2}\), and thus \(\displaystyle \theta \left( x\le \displaystyle \frac{y}{2} \right) =1\). In the first step, we change the order of integration

$$\begin{aligned} \int _0^\infty v\mathscr {B}_2(v){\textrm{d}}v&= \int _0^\infty v \int _v^\infty n(t, y)\int _0^{y} \Omega _t\left( y\right) \beta \left( x, y\right) \delta \left[ x- \left( y-v\right) \right] \\&\qquad \theta \left( x\le \displaystyle \frac{y}{2}\right) {\textrm{d}}x{\textrm{d}}y {\textrm{d}}v \\&= \int _0^\infty n(t, y)\Omega _t\left( y\right) \int _0^y \int _0^{\displaystyle \frac{y}{2}} v\beta (x,y)\delta \left[ x-\left( y-v\right) \right] {\textrm{d}}x{\textrm{d}}v {\textrm{d}}y \\&= \int _0^\infty n(t, y)\Omega _t\left( y\right) \int _0^y \int _0^{\displaystyle \frac{y}{2}} v\beta (x,y)\delta \left[ v-\left( y-x\right) \right] {\textrm{d}}x{\textrm{d}}v {\textrm{d}}y \\&= \int _0^\infty n(t, y)\Omega _t\left( y\right) \int _{\displaystyle \frac{y}{2}}^y v\beta (v,y) {\textrm{d}}v {\textrm{d}}y. \end{aligned}$$

Therefore,

$$\begin{aligned} \int _0^\infty v\mathscr {B}(v) {\textrm{d}}v&= \int _0^\infty v\mathscr {B}_1(v){\textrm{d}}v + \int _0^\infty v\mathscr {B}_2(v){\textrm{d}}v \\&= \int _0^\infty n(t, y)\Omega _t(y) \left[ \int _0^{\displaystyle \frac{y}{2}} v\beta (v,y) {\textrm{d}}v\right] {\textrm{d}}y + \int _0^\infty n(t, y)\Omega _t\left( y\right) \\&\quad \int _{\displaystyle \frac{y}{2}}^y v\beta (v,y) {\textrm{d}}v {\textrm{d}}y \\&= \int _0^\infty n(t, y)\Omega _t(y) \int _0^y v\beta (v,y) {\textrm{d}}v {\textrm{d}}y \\&= \int _0^\infty yn(t, y)\Omega _t(y) {\textrm{d}}y. \end{aligned}$$

Also

$$\begin{aligned} \int _0^\infty v\mathscr {D}(v) {\textrm{d}}v = \int _0^\infty yn(t, y)\Omega _t(y) {\textrm{d}}y. \end{aligned}$$

Therefore the mass conservation law follows. Similarly we can calculate the zeroth moment as

$$\begin{aligned} \frac{{\textrm{d}}M_0(t)}{{\textrm{d}}t} = \int _0^\infty \left[ \mathscr {B}(v) - \mathscr {D}(v)\right] {\textrm{d}}v = \int _0^\infty n(t, y)\Omega _t(y) {\textrm{d}}y \end{aligned}$$

1.3 State-of-the-art and motivation

Due to the complex nature of the fragmentation equation, few analytical solutions are available even for simple fragmentation kernels [2, 28,29,30]. For complex kernels, various numerical methods have been proposed during the last few decades including the fixed pivot technique [8], the cell average technique [9] and finite volume schemes [7, 17, 21, 23, 25]. Some other numerical methods related to multidimensional fragmentation equations can be found in [16, 18, 20, 22]. The merits and shortcomings of the fixed pivot technique and cell average methods have been discussed in detail by [9, 24]. It was shown that the fixed pivot technique is second order convergent on uniform and non-uniform grids whereas first order convergence on locally uniform grids. Moreover, for non-uniform random grids, the fixed pivot technique is inconsistent. In addition, Kumar and Warnecke [9] demonstrated that the cell average technique is second order convergent on uniform, non-uniform and locally uniform grids. The cell average technique was shown to exhibit first order convergence on non-uniform random grid and perform better than the fixed pivot technique. The main disadvantage of these two techniques is that both methods are inefficient in the case of partial breakup kernels [12, 14].

Furthermore, it was demonstrated that the recently developed finite volume methods are more efficient than both the fixed pivot and the cell average technique, however they require modification of the mathematical formulations to accommodate partial breakup kernels. However, regardless of the choice of grid existing finite volume schemes exhibit second order convergence. The recently developed fixed pivot technique [11] overcomes all of the shortcomings of the aforementioned numerical methods and seems to be more flexible in terms of implementing any kind of breakup kernel. Furthermore, because the breakup rate in a CFD simulation is dependent on flow field parameters, the integration (1.4) would have to be performed for each control volume, and that significantly increases the computational cost [11]. Due to these reasons, the modified fixed pivot technique [11] is more powerful and flexible for tackling real life applications in chemical engineering and pharmaceutical sciences. However, the convergence analysis of the recently developed fixed pivot technique for a binary fragmentation equation is still incomplete and continues to be a major problem in the literature. The complexity in conducting the convergence analysis of the recently developed fixed pivot technique is due to the presence of the two tedious factors (2.1) and (2.2) related to the distribution of the particle properties to the neighbouring nodes. In order to show the reliability of the method, the comparison with the existing cell average technique [9] is conducted. Moreover, the experimental order of convergence is used to confirm the theoretical rate of convergence for different grids.

The rest of the current article is structured as follows: Sect. 2 is devoted for providing the mathematical formulation of the fixed pivot technique along with the theoretical proofs of mass conservation law and number preservation. Next Sect. 3 focuses on the convergence and stability analysis of the fixed pivot technique while discussing consistency of the method in detail. Further, numerical results in terms of moments, number of particles, average size particles and experimental order of convergence are compared and discussed in Sect. 4. Finally, some remarks, conclusions and future aspects are discussed in Sect. 5.

2 Computational domain and numerical scheme

Let the computational domain \(\displaystyle \left[ 0, x_{\max }\right] \) be subdivided into I subintervals \(\displaystyle \Lambda _i := \left[ b_i, b_{i+1} \right] .\) The midpoint of each cell \(\Lambda _i\) is denoted by \(x_i\), and the cell width is denoted by \(\displaystyle \Delta x_i := b_{i+1}-b_i\). Accordingly, the discrete number of particles in each cell (\(N_i\)) calculated from the number density function reads as

$$\begin{aligned} N_i \approx \int _{b_i}^{b_{i+1}} n(t, v){\textrm{d}}v, \end{aligned}$$

To derive the mathematical formulation of the scheme, the birth term is integrated over the size range \(\Lambda _i\),

$$\begin{aligned} \int _{b_i}^{b_{i+1}} \mathscr {B}(v){\textrm{d}}v =&\int _{b_i}^{b_{i+1}} \int _v^\infty n(t, y)\Omega _t\left( y\right) \beta (x, y)\theta \left( v\le \frac{y}{2}\right) {\textrm{d}}y {\textrm{d}}v \\&+ \int _{b_i}^{b_{i+1}} \int _v^\infty n(t, y) \int _0^{y} \Omega _t\left( y\right) \beta (x, y)\delta \left[ x- \left( y-v\right) \right] \\&\quad \theta \left( x\le \frac{y}{2}\right) {\textrm{d}}x{\textrm{d}}y {\textrm{d}}v. \end{aligned}$$

Applying several approximations, the birth term of the modified fixed pivot technique is written as

$$\begin{aligned} \hat{\mathscr {B}}\left( x_i\right) =&\sum _{j=i}^I \hat{N}_j\Omega _t\left( x_j\right) \beta \left( x_i,x_j\right) \int _{b_i}^{b_{i+1}} \theta \left( v\le \frac{x_j}{2}\right) {\textrm{d}}v \\&+ \sum _{j=i}^I \sum _{k=1}^j \hat{N}_j\Omega _t\left( x_j\right) \beta \left( x_k, x_j\right) \int _{b_i}^{b_{i+1}} \delta \left[ x_k - \left( x_j-v\right) \right] {\textrm{d}}v \\&\quad \int _{b_k}^{b_{k+1}} \theta \left( x\le \frac{x_j}{2} \right) {\textrm{d}}x \\ =&\sum _{j=i}^I \hat{N}_j\left[ \Omega _t\left( x_j\right) \beta \left( x_i, x_j\right) \Delta \nu _i(j) + \sum _{k=1}^j \Omega _t\left( x_j\right) \beta \left( x_k, x_j\right) Y_{ijk}\Delta \nu _k(j)\right] \\ =:\,&\hat{\mathscr {B}}_1\left( x_i\right) + \hat{\mathscr {B}}_2\left( x_i\right) . \end{aligned}$$

Here, \(\hat{N}_i\) is the numerical approximation of the number of particles in each cell, \(\Delta \nu _i(j)\), \(Y_{ijk}\) are the birth-modification factors which assign the particles to their nearest representatives, and are defined as

$$\begin{aligned} \Delta \nu _i(j)&:= \left\{ \begin{array}{lll} \displaystyle b_{i+1} - b_i, &{} \text{ when } \displaystyle b_{i+1}\le \frac{x_j}{2}, \\ \displaystyle \frac{x_j}{2} - b_i, &{} \text{ when } \displaystyle b_i \le \frac{x_j}{2} < b_{i+1}, \\ \displaystyle 0, &{}\text{ when } \displaystyle b_i \ge \frac{x_j}{2} , \end{array} \right. \end{aligned}$$
(2.1)
$$\begin{aligned} Y_{ijk}&:= \left\{ \begin{array}{lll} \displaystyle \frac{\left( x_j-x_k\right) -x_{i-1}}{x_i - x_{i-1}}, &{} \text{ when } \displaystyle x_{i-1}\le x_j-x_k<x_i, \\ \displaystyle \frac{x_{i+1}-\left( x_j-x_k\right) }{x_{i+1}-x_i}, &{} \text{ when } \displaystyle x_i \le x_j - x_k < x_{i+1}, \\ \displaystyle 0, &{}\text{ otherwise. } \end{array} \right. \end{aligned}$$
(2.2)

Similarly, integrating the death term over each subinterval gives

$$\begin{aligned} \hat{\mathscr {D}}\left( x_i\right) := \hat{N}_i\sum _{j=1}^i \Omega _t\left( x_j\right) \beta \left( x_i, x_j\right) \Delta \nu _i(j). \end{aligned}$$

Note that Liao et al. [11] have written the discrete partial breakage kernel \(\Omega _p\) as follows:

$$\begin{aligned} \Omega _p\left( x_i, x_j\right) := \Omega _t\left( x_j\right) \beta \left( x_i, x_j\right) \quad \text{ and } \quad \Omega _t\left( x_i\right) := \sum _{j=1}^i \Omega _p\left( x_j, x_i\right) \Delta \nu _i(j). \end{aligned}$$

Under the above consideration, the modified fixed pivot technique of [11] is written as

$$\begin{aligned} \frac{{\textrm{d}}\hat{N}_i}{{\textrm{d}}t}&= \underbrace{\sum _{j=i}^I \hat{N}_j\left[ \Omega _p\left( x_i, x_j\right) \Delta \nu _i(j) + \sum _{k=1}^j \Omega _p\left( x_k, x_j\right) Y_{ijk}\Delta \nu _{k}(j) \right] }_{= \hat{\mathscr {B}}\left( x_i\right) }\nonumber \\&\quad - \underbrace{\hat{N}_i\sum _{j=1}^i \Omega _p\left( x_j, x_i\right) \Delta \nu _j(i)}_{= \hat{\mathscr {D}}\left( x_i\right) }. \end{aligned}$$
(2.3)

2.1 Zeroth order moment consistency

During fragmentation of \(x_j\), the distribution of daughter particles is based upon their sizes relative to \(\displaystyle \frac{x_j}{2}\). Consider that \(\displaystyle \frac{x_j}{2}\) falls in some cell \(\displaystyle \Lambda _{j'}\) where \(1\le j'\le j\). Then \(\hat{\mathscr {B}}_1\left( x_i\right) \) and \(\hat{\mathscr {B}}_2\left( x_i\right) \), respectively, gives the distribution of particles smaller and larger than the half-range size of \(x_j\). In this regard, we first consider the particle size distribution smaller than \(\displaystyle \frac{x_j}{2}\) and get

$$\begin{aligned} \sum _{i=1}^I \hat{\mathscr {B}}_1\left( x_i\right)&= \sum _{i=1}^I \sum _{j=i}^I \hat{N}_j\Omega _t\left( x_j\right) \beta \left( x_i,x_j\right) \int _{b_i}^{b_{i+1}} \theta \left( v\le \frac{x_j}{2}\right) {\textrm{d}}v \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^j \beta \left( x_i,x_j\right) \int _{b_i}^{b_{i+1}} \theta \left( v\le \frac{x_j}{2}\right) {\textrm{d}}v. \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^{j'} \beta \left( x_i,x_j\right) \Delta \nu _i(j). \end{aligned}$$

The distribution of particles larger than \(\displaystyle \frac{x_j}{2}\) is achived through the complementary particles \(\displaystyle x \left( \le \frac{x_j}{2}\right) \). Due to symmetry of the breakage function, x falls in the same interval \(1\le j'\le j\), and thus we get

$$\begin{aligned} \sum _{i=1}^I \hat{\mathscr {B}}_2\left( x_i\right)&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^j \sum _{k=1}^j \beta \left( x_k,x_j\right) \int _{b_i}^{b_{i+1}} \delta \left[ x_k-\left( x_j-v\right) \right] {\textrm{d}}v\\&\quad \int _{b_k}^{b_{k+1}} \theta \left( x\le \frac{x_j}{2}\right) {\textrm{d}}x. \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^{j} \sum _{k=1}^{j'} \beta \left( x_k,x_j\right) \Delta \nu _k(j) \int _{b_i}^{b_{i+1}} \delta \left[ v-\left( x_j-x_k\right) \right] {\textrm{d}}v\\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=j'}^{j} \beta \left( x_i,x_j\right) \Delta \nu _i(j). \end{aligned}$$

Combining the two cases, we get

$$\begin{aligned} \sum _{i=1}^I \hat{\mathscr {B}}\left( x_i\right)&= \sum _{i=1}^I \hat{\mathscr {B}}_1\left( x_i\right) + \sum _{i=1}^I \hat{\mathscr {B}}_2\left( x_i\right) = \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^{j} \beta \left( x_i,x_j\right) \Delta x_i \\&= 2\sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) . \end{aligned}$$

Hence, the time evolution of the discrete zeroth moment is calculated as

$$\begin{aligned} \frac{d\hat{M}_0}{dt} = \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) . \end{aligned}$$

2.2 Mass conservation law

First consider the distribution of the particles of size smaller than \(\displaystyle \frac{x_j}{2}\).

$$\begin{aligned} \sum _{i=1}^I x_i \hat{\mathscr {B}}_1\left( x_i\right)&= \sum _{i=1}^I x_i \sum _{j=i}^I \hat{N}_j\Omega _t\left( x_j\right) \beta \left( x_i,x_j\right) \int _{b_i}^{b_{i+1}} \theta \left( v\le \frac{x_j}{2}\right) {\textrm{d}}v \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^j x_i\beta \left( x_i,x_j\right) \int _{b_i}^{b_{i+1}} \theta \left( v\le \frac{x_j}{2}\right) {\textrm{d}}v. \end{aligned}$$

Now the argument of cell allocation for particles follows similarly to the previous description, that is, consider that \(\displaystyle \frac{x_j}{2}\) falls in some cell \(\displaystyle \Lambda _{j'}\) such that \(1\le j'\le j\). Therefore,

$$\begin{aligned} \sum _{i=1}^I x_i \hat{\mathscr {B}}_1\left( x_i\right) = \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^{j'} x_i \beta \left( x_i,x_j\right) \Delta \nu _i(j). \end{aligned}$$

For the distribution of particles larger than \(\displaystyle \frac{x_j}{2}\), we compute

$$\begin{aligned} \sum _{i=1}^I x_i \hat{\mathscr {B}}_2\left( x_i\right)&= \sum _{i=1}^I x_i \sum _{j=i}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{k=1}^j \beta \left( x_k,x_j\right) \int _{b_i}^{b_{i+1}} \delta \left[ x_k-\left( x_j-v\right) \right] {\textrm{d}}v\\&\quad \int _{b_k}^{b_{k+1}} \theta \left( x\le \frac{x_j}{2}\right) {\textrm{d}}x \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^j x_i \sum _{k=1}^j \beta \left( x_k,x_j\right) \int _{b_i}^{b_{i+1}} \delta \left[ x_k-\left( x_j-v\right) \right] {\textrm{d}}v\\&\quad \times \int _{b_k}^{b_{k+1}} \theta \left( x\le \frac{x_j}{2}\right) {\textrm{d}}x. \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^{j} x_i \sum _{k=1}^{j'} \beta \left( x_k,x_j\right) \Delta \nu _k(j)\\&\quad \times \int _{b_i}^{b_{i+1}} \delta \left[ v-\left( x_j-x_k\right) \right] {\textrm{d}}v \\&= \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=j'}^{j} x_i \beta \left( x_i,x_j\right) \Delta \nu _i(j). \end{aligned}$$

Combining the two cases, we get

$$\begin{aligned} \sum _{i=1}^I x_i \hat{\mathscr {B}}\left( x_i\right)&= \sum _{i=1}^I x_i \hat{\mathscr {B}}_1\left( x_i\right) + \sum _{i=1}^I x_i \hat{\mathscr {B}}_2\left( x_i\right) = \sum _{j=1}^I \hat{N}_j\Omega _t\left( x_j\right) \sum _{i=1}^{j} x_i \beta \left( x_i,x_j\right) \!\Delta x_i \\&= \sum _{j=1}^I x_j \hat{N}_j\Omega _t\left( x_j\right) . \end{aligned}$$

Hence, time evolution of the discrete first moment is calculated as

$$\begin{aligned} \frac{d\hat{M}_1}{dt} = 0. \end{aligned}$$

3 Convergence analysis

Consider the numerical solution expressed in the vector form \(\displaystyle \hat{\textbf{N}} = \{\hat{N}_1, \hat{N}_2, \ldots , \hat{N}_I \}\). Then the semi-discrete system (2.3) is expressed as

$$\begin{aligned} \frac{{\textrm{d}}\hat{\textbf{N}}}{{\textrm{d}}t} = \mathscr {A}\hat{\textbf{N}}, \end{aligned}$$
(3.1)

where \(\mathscr {A}\) is the \(I\times I\) upper triangular matrix

$$\begin{aligned} \mathscr {A} := \begin{bmatrix} \alpha _{11} - \Omega _t\left( x_1\right) &{} \alpha _{12} &{} \alpha _{13} &{} \cdots &{} \alpha _{1I} \\ 0 &{} \alpha _{22} - \Omega _t\left( x_2\right) &{} \alpha _{23} &{} \cdots &{} \alpha _{2I} \\ \vdots &{} \vdots &{} \vdots &{} \ddots &{} \vdots \\ 0 &{} 0 &{} 0 &{} \cdots &{} \alpha _{II} - \Omega _t\left( x_I\right) \end{bmatrix}, \end{aligned}$$

with the components \( \displaystyle \alpha _{ij} := \Omega _t\left( x_j\right) \beta \left( x_i,x_j\right) \Delta \nu _i\,(j) + \sum \nolimits _{k=1}^j \Omega _t\left( x_j\right) \beta \left( x_k,x_j\right) Y_{ijk} \Delta \nu _k(j)\).

Throughout the subsequent discussion, we work with the discrete \(L^1-\)norm defined as follows

$$\begin{aligned} \left\| \hat{N}(t)\right\| = \sum _{i=1}^I \left| \hat{N}_i(t)\right| . \end{aligned}$$

The consistency and stability of the model are analyzed with the help of following definitions and theorems which are quoted from Hundsdorfer and Verwer [4].

Definition 3.1

The spatial truncation error is defined by the residual resulting from substituting the exact solution \(\displaystyle \textbf{N}= \{N_1, N_2, \ldots , N_I \}\) into the discrete system as

$$\begin{aligned} \sigma (t) = \frac{{\textrm{d}}\textbf{N}(t)}{{\textrm{d}}t} - \mathscr {A}\textbf{N}(t). \end{aligned}$$

The scheme is called consistent of order q if, for \(\Delta x \rightarrow 0\),

$$\begin{aligned} \left\| \sigma (t) \right\| = \mathcal {O}\left( \Delta x^q\right) \quad \text{ uniformly } \text{ for } \text{ all }\quad 0\le t\le T. \end{aligned}$$

Definition 3.2

Consider a matrix \(\displaystyle A := \left[ a_{ij}\right] \in {\mathbb R}^{m\times m}\). Then for the \(L^p\) vector norm, the corresponding logarithmic norm of A is defined as

$$\begin{aligned} \tilde{\mu }(A) = \lim _{\tau \rightarrow 0} \frac{\left\| I+\tau A\right\| -1}{\tau }. \end{aligned}$$

Moreover, the logarithmic norm of a matrix A corresponding to the \(L^1\) norm on \({\mathbb R}^m\) is given by

$$\begin{aligned} \tilde{\mu }_1(A) = \max _j \left[ a_{jj} + \sum _{i\ne j}\left| a_{ij}\right| \right] . \end{aligned}$$

Theorem 3.1

If \(A := \left[ a_{ij}\right] \in {\mathbb R}^{m\times m}\) and \(\omega \in {\mathbb R}\) then

$$\begin{aligned} \tilde{\mu }(A) \le \omega \quad \text{ if } \text{ and } \text{ only } \text{ if }\quad \left\| \exp (tA)\right\| \le \left\| \exp (t\omega )\right\| \quad \text{ for } \text{ all } \quad t \ge 0. \end{aligned}$$

Definition 3.3

The semi-discrete system \(\displaystyle \frac{{\textrm{d}}w(t)}{{\textrm{d}}t} = Aw(t)\) is stable if

$$\begin{aligned} \left\| \exp (tA)\right\| \le K\exp \left( t\omega \right) \quad \text{ for } \text{ all } \quad 0\le t \le T \end{aligned}$$

holds for all grids with constant \(K\ge 1\) and \(\omega \in {\mathbb R}\), both independent of \(\Delta x\).

Theorem 3.2

The solution of the linear semi-discrete system \(\displaystyle \frac{{\textrm{d}}w(t)}{{\textrm{d}}t} = Aw(t)\) is nonnegative if and only if

$$\begin{aligned} a_{ij} \ge 0, \quad \text{ for } \text{ all } \quad i\ne j, \end{aligned}$$

where \(a_{ij}\) are the real entries of the matrix A.

Proof

The proof of the above theorem can be found in [4] (Theorem 7.2, p. 117). \(\square \)

3.1 Nonnegativity

Lemma 3.3

(Nonnegativity) The numerical solution of the semi-discrete scheme (3.1) is nonnegative.

Proof

From the definitions of all the functions and the matrix \(\mathscr {A}\), it can be easily seen that \(\alpha _{ij}\ge 0\) for all \(i\ne j\). Thus recalling Theorem 3.2, we deduce that the solution of the semi-discrete system (3.1) is nonnegative. \(\square \)

3.2 Consistency

Let \(\mathcal {C}^2\left[ a,b\right] \) denote the space of twice continuously differentiable functions on the nonempty bounded closed interval [a, b].

Lemma 3.4

(Consistency) Suppose that \(\displaystyle \beta \in \mathcal {C}^2\left( \left[ 0, x_{\max }\right] \times \left[ 0, x_{\max }\right] \right) \) and \(\displaystyle \Omega _t\in \mathcal {C}^2\left( \left[ 0, x_{\max }\right] \right) \). The semi-discrete scheme (3.1) is second order accurate, independent of the choice of the spatial mesh used.

Proof

The spatial truncation error is given by

$$\begin{aligned} \sigma _i(t) = \frac{{\textrm{d}}\textbf{N}}{{\textrm{d}}t} - \frac{{\textrm{d}}\hat{\textbf{N}}}{{\textrm{d}}t} = \mathscr {A}\textbf{N}- \mathscr {A}\hat{\textbf{N}} = \left[ \mathscr {B}(x_i) - \hat{\mathscr {B}}(x_i)\right] - \left[ \mathscr {D}\left( x_i\right) - \hat{\mathscr {D}}\left( x_i\right) \right] . \end{aligned}$$

Integrating the birth term \(\displaystyle \mathscr {B}\left( v\right) \) over each cell \(\Lambda _i\) we get that

$$\begin{aligned} \int _{b_i}^{b_{i+1}} \mathscr {B}\left( v\right) {\textrm{d}}v =&\int _{b_i}^{b_{i+1}} \int _v^{x_{\max }} n(t, y)\Omega _t\left( y\right) \beta (x, y)\theta \left( v\le \frac{y}{2}\right) {\textrm{d}}y {\textrm{d}}v \\&+ \int _{b_i}^{b_{i+1}} \int _v^{x_{\max }} n(t, y) \int _0^{y} \Omega _t\left( y\right) \beta (x, y)\delta \left[ x- \left( y-v\right) \right] \\&\quad \theta \left( x\le \frac{y}{2}\right) {\textrm{d}}x{\textrm{d}}y {\textrm{d}}v. \end{aligned}$$

Applying quadrature rules for the outer integrals on the right-hand side, we get

$$\begin{aligned} \mathscr {B}\left( x_i\right) =&\int _{x_i}^{x_{\max }} n(t, y)\Omega _t\left( y\right) \beta (x_i, y)\theta \left( x_i\le \frac{y}{2}\right) \Delta x_i {\textrm{d}}y \\&+ \int _{x_i}^{x_{\max }} n(t, y) \int _0^{y} \Omega _t\left( y\right) \beta (x, y)\delta \left[ x- \left( y-x_i\right) \right] \theta \left( x\le \frac{y}{2}\right) \Delta x_i {\textrm{d}}x{\textrm{d}}y + \mathcal {O}\left( \Delta x_i^3\right) \\ =&\sum _{k=i}^I \int _{\Lambda _k} n(t, y)\Omega _t\left( y\right) \beta (x_i, y)\theta \left( x_i\le \frac{y}{2}\right) \Delta x_i {\textrm{d}}y \\&+ \sum _{k=i}^I \int _{\Lambda _k} n(t, y) \int _0^{y} \Omega _t\left( y\right) \beta (x, y)\delta \left[ x- \left( y-x_i\right) \right] \theta \left( x\le \frac{y}{2}\right) \Delta x_i {\textrm{d}}x{\textrm{d}}y \\&- \underbrace{\int _{b_i}^{x_i} n(t, y)\Omega _t\left( y\right) \beta (x_i, y)\theta \left( x_i\le \frac{y}{2}\right) \Delta x_i {\textrm{d}}y}_{E_1} \\&- \underbrace{\int _{b_i}^{x_i} n(t, y)\int _0^{y} \Omega _t\left( y\right) \beta (x, y)\delta \left[ x- \left( y-x_i\right) \right] \theta \left( x\le \frac{y}{2}\right) \Delta x_i {\textrm{d}}x{\textrm{d}}y}_{E_2} + \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

Again we apply the midpoint rule to the outer integrals of the first two terms and get

$$\begin{aligned} \mathscr {B}\left( x_i\right) =&\sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \beta (x_i, x_k)\theta \left( x_i\le \frac{x_k}{2}\right) \Delta x_i \\&+ \sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \int _0^{x_k} \beta (x, x_k)\delta \left[ x- \left( x_k-x_i\right) \right] \theta \left( x\le \frac{x_k}{2}\right) \Delta x_i {\textrm{d}}x\\&- E_1 - E_2 + \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

Before we proceed further, let us first gather the estimates of a few terms which will be required in the subsequent discussions. Application of midpoint quadrature rules gives

$$\begin{aligned} \int _{\Lambda _i} \theta \left( x\le \frac{x_k}{2}\right) {\textrm{d}}x = \theta \left( x_i\le \frac{x_k}{2}\right) \Delta x_i + \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

Therefore,

$$\begin{aligned} \mathscr {B}\left( x_i\right) =&\sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \beta (x_i, x_k)\int _{\Lambda _i} \theta \left( x\le \frac{x_k}{2}\right) {\textrm{d}}x \\&+ \sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \sum _{j=1}^k \int _{\Lambda _j} \beta (x, x_k)\delta \left[ x- \left( x_k-x_i\right) \right] \theta \left( x\le \frac{x_k}{2}\right) \Delta x_i {\textrm{d}}x \\&- E_1 - E_2 - E_3 + \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

Here,

$$\begin{aligned} E_3 = \sum _{k=i}^I \hat{N}_k\Omega _t\left( x_k\right) \int _{x_k}^{b_{k+1}} \beta (x, x_k)\delta \left[ x- \left( x_k-x_i\right) \right] \theta \left( x\le \frac{x_k}{2}\right) \Delta x_i {\textrm{d}}x. \end{aligned}$$

By the application of the midpoint rule, the definition of \(\displaystyle \Delta \nu _i(j)\), \(Y_{ijk}\) and other simple computations, we get

$$\begin{aligned} \mathscr {B}\left( x_i\right) =&\sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \beta (x_i, x_k)\int _{\Lambda _i} \theta \left( x\le \frac{x_k}{2}\right) {\textrm{d}}x \\&+ \sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \sum _{j=1}^k \beta (x_j, x_k)\! \int _{\Lambda _i} \!\delta \left[ x_j - \left( x_k- v\right) \right] {\textrm{d}}v \int _{\Lambda _j} \theta \left( x\le \frac{x_k}{2}\right) {\textrm{d}}x \\&- E_1 - E_2 - E_3 + \mathcal {O}\left( \Delta x_i^3\right) \\ =&\sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \beta (x_i, x_k)\Delta \nu _i(k) + \sum _{k=i}^I \hat{N}_k \Omega _t\left( x_k\right) \sum _{j=1}^k \beta (x_j, x_k) Y_{ijk} \Delta \nu _j(k) \\&- E_1 - E_2 - E_3 + \mathcal {O}\left( \Delta x_i^3\right) \\ = \,&\hat{\mathscr {B}}\left( x_i\right) - E_1 - E_2 - E_3 + \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

We now proceed to estimate the error terms \(E_1\), \(E_2\) and \(E_3\). In general, the definition of \(\beta \) yields that

$$\begin{aligned} \beta (x,y) = 0 \quad \text{ whenever } \quad x\ge y. \end{aligned}$$

Therefore, the definition of \(E_1\) gives

$$\begin{aligned} \int _{b_i}^{x_i} \beta \left( x_i, y\right) {\textrm{d}}y = 0 \quad \text{ which } \text{ implies } \quad E_1 =0. \end{aligned}$$

Similarly,

$$\begin{aligned} \int _{x_k}^{b_{k+1}} \beta \left( x, x_k\right) {\textrm{d}}x = 0 \quad \text{ which } \text{ implies } \quad E_3 =0. \end{aligned}$$

Now for estimating \(E_2\), limits of the integrals give

$$\begin{aligned} b_i \le y \le x_i \implies b_i - x_i \le y - x_i \le 0 \quad \text{ and } \quad 0\le x \le y, \end{aligned}$$

which then imply that

$$\begin{aligned} \delta \left[ x- \left( y-x_i\right) \right] = 0 \quad \text{ that } \text{ is } \quad E_2 =0. \end{aligned}$$

Hence, we get

$$\begin{aligned} \mathscr {B}\left( x_i\right) - \hat{\mathscr {B}}\left( x_i\right) = \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

For the death term, direct application of midpoint quadrature formula leads to

$$\begin{aligned} \mathscr {D}\left( x_i\right) - \hat{\mathscr {D}}\left( x_i\right) = \mathcal {O}\left( \Delta x_i^3\right) . \end{aligned}$$

Hence,

$$\begin{aligned} \sigma _i(t) = \mathcal {O}\left( \Delta x_i^3\right) \quad \text{ which } \text{ implies } \quad \Vert \sigma (t)\Vert = \mathcal {O}\left( \Delta x^2\right) . \end{aligned}$$

This proves that the scheme is second order consistent irrespective of the choice of the spatial mesh. \(\square \)

3.3 Stability

Proof

In order to establish the stability of the scheme, we now compute the logarithmic norm (Def. 3.2) of the matrix A as

$$\begin{aligned} \tilde{\mu }_1(A) = \max _j\left[ a_{jj} + \sum _{i\ne j} |a_{ij}| \right] . \end{aligned}$$

Since all non-diagonal entries of the matrix \(\mathscr {A}\) are nonnegative, the above logarithmic norm takes the following form

$$\begin{aligned} \tilde{\mu }_1(\mathscr {A}) = \max _j\left[ \sum _{i} |a_{ij}|\right] = \max _j\left[ \sum _{i=1}^I |a_{ij}|\right] = \max _j\left[ \sum _{i=1}^j \alpha _{ij} - \Omega _t\left( x_j\right) \right] . \end{aligned}$$

Now,

$$\begin{aligned} \sum _{i=1}^j \alpha _{ij}&= \sum _{i=1}^j \Omega _t\left( x_j\right) \beta \left( x_i,x_j\right) \Delta \nu _i(j) + \sum _{i=1}^j\sum _{k=1}^j \Omega _t\left( x_j\right) \beta \left( x_k,x_j\right) Y_{ijk} \Delta \nu _k(j) \\&= \Omega _t\left( x_j\right) \sum _{i=1}^j \beta \left( x_i,x_j\right) \Delta \nu _i(j) + \Omega _t\left( x_j\right) \sum _{i=1}^j \sum _{k=1}^j \beta \left( x_k,x_j\right) Y_{ijk} \Delta \nu _k(j). \end{aligned}$$

Suppose that \(\displaystyle \frac{x_j}{2}\) falls in some interval \(\Lambda _{j'}\) satisfying \(1\le j' \le j\). Then the definition of \(\Delta \nu _i(j)\) redefines the formulation as

$$\begin{aligned} \sum _{i=1}^j \alpha _{ij}&= \Omega _t\left( x_j\right) \sum _{i=1}^{j'} \beta \left( x_i,x_j\right) \Delta \nu _i(j) + \Omega _t\left( x_j\right) \sum _{i=1}^j \sum _{k=1}^{j'} \beta \left( x_k,x_j\right) Y_{ijk} \Delta \nu _k(j). \end{aligned}$$

Again from the definition, the lower boundary of \(Y_{ijk}\) is given by \(x_{i-1}\le x_j-x_k\le x_i\) which implies that the indices i and k are together replaced by a new index \(i'\) which ranges from \(j'+1 \rightarrow j\). Therefore, further replacing \(i'\) by i, we get

$$\begin{aligned} \sum _{i=1}^j \alpha _{ij}&= \Omega _t\left( x_j\right) \sum _{i=1}^{j'} \beta \left( x_i,x_j\right) \Delta \nu _i(j) + \Omega _t\left( x_j\right) \sum _{i=j'+1}^j \beta \left( x_j-x_i,x_j\right) \Delta \nu _i(j) \\&\le \Omega _t\left( x_j\right) \sum _{i=1}^{j} \beta \left( x_i,x_j\right) \Delta x_i \\&\le 2\max _{x\in [0,x_{\max }]} \Omega _t(x) = 2\omega \quad \text{(say) }. \end{aligned}$$

Therefore the logarithmic norm \(\tilde{\mu }_1(A)\le 2\omega \), and consequently by using Theorem 3.1 we get \(\Vert \exp (tA)\Vert \le \exp (2\omega t)\) which ensures the stability of the scheme with stability constant \(2\omega \). \(\square \)

Lipschitz condition:

Consider the problem written as

$$\begin{aligned} \frac{{\textrm{d}}\hat{\textbf{N}}}{{\textrm{d}}t} = \textbf{J}(\hat{\textbf{N}}), \end{aligned}$$

where the numerical flux is defined by

$$\begin{aligned} \textbf{J}(\hat{\textbf{N}})&:= \sum _{j=i}^I \hat{N}_j\left[ \Omega _t\left( x_j\right) \beta \left( x_i, x_j\right) \Delta \nu _i(j) + \sum _{k=1}^j \Omega _t\left( x_j\right) \beta \left( x_k, x_j\right) Y_{ijk}\Delta \nu _k(j)\right] \\&\quad - \hat{N}_i\Omega _t\left( x_i\right) . \end{aligned}$$

Theorem 3.5

(Ref. [13]) If \(\textbf{J}(\textbf{N})\) satisfies the Lipschitz condition

$$\begin{aligned} \left\| \textbf{J}(\textbf{N}) - \textbf{J}(\hat{\textbf{N}}) \right\| \le \xi \left\| \textbf{N}- \hat{\textbf{N}}\right\| ,\quad \xi <\infty , \end{aligned}$$

for all \(0\le t \le T\) and \(\textbf{N},\hat{\textbf{N}}\in {\mathbb R}^I\), then a consistent discretization is also convergent. Furthermore, the order of convergence is the same as the order of consistency.

Therefore,

Changing the order of the summation, we get

$$\begin{aligned} \left\| \textbf{J}(\textbf{N}) - \textbf{J}(\hat{\textbf{N}}) \right\| \le&\sum _{j=1}^I \left| N_j - \hat{N}_j \right| \Omega _t\left( x_j\right) \left[ \sum _{i=1}^j \beta \left( x_i,x_j\right) \Delta \nu _i(j) \right. \\&\left. + \sum _{i=1}^j \sum _{k=1}^{j} \beta \left( x_k,x_j\right) Y_{ijk} \Delta _k(j) \right] \\&+ \sum _{i=1}^I \left| N_i - \hat{N}_i \right| \Omega _t\left( x_i\right) . \end{aligned}$$

Supposing that \(\displaystyle \frac{x_j}{2}\) falls on some interval \(1\le j'\le j\), we get

$$\begin{aligned} \left\| \textbf{J}(\textbf{N}) - \textbf{J}(\hat{\textbf{N}}) \right\| \le&\sum _{j=1}^I \left| N_j - \hat{N}_j \right| \Omega _t\left( x_j\right) \left[ \sum _{i=1}^{j'} \beta \left( x_i,x_j\right) \Delta \nu _i(j) \right. \\&\left. + \sum _{i=1}^j \sum _{k=1}^{j'} \beta \left( x_k,x_j\right) Y_{ijk} \Delta _k(j) \right] + \sum _{i=1}^I \left| N_i - \hat{N}_i \right| \Omega _t\left( x_i\right) \\ \le&\sum _{j=1}^I \left| N_j - \hat{N}_j \right| \Omega _t\left( x_j\right) \left[ \sum _{i=1}^{j'} \beta \left( x_i,x_j\right) \Delta x_i \right. \\&\left. + \sum _{i=j'+1}^{j} \beta \left( x_i,x_j\right) \Delta x_i \right] + \sum _{i=1}^I \left| N_i - \hat{N}_i \right| \Omega _t\left( x_i\right) \\ \le&3\omega \left\| \textbf{N}- \hat{\textbf{N}} \right\| , \end{aligned}$$

where \(\displaystyle \omega := \max _{x\in [0,x_{\max }]} \Omega _t(x)\). Therefore the Lipschitz constant \(\xi := 3\omega <\infty \). \(\square \)

4 Numerical results and discussion

This part of the paper is devoted to checking the accuracy of the modified fixed pivot technique (MFPT) against the cell average technique (CAT) [9] and exact results [29, 30] for different initial conditions. The detailed comparison of the traditional fixed pivot technique [10] against the modified one is provided in Liao et al. [11] for different selection functions and binary breakage kernels. It has been shown that the MFPT performs much better than the traditional fixed pivot technique [10]. In addition, Liao et al. [11] have also shown that the MFPT can be easily adjusted to incorporate partial breakup kernels which is one of the main requirements when solving problems in chemical and pharmaceutical sciences using CFD tools. So we omit these results in our study.

Due to the extensive usage of the CAT for solving real life applications such as depolymerization [1] and granulation [5], the comparison of the MFPT is extended against the CAT. The optimization of various breakage kernel and selection function parameters is essential for resolving the accurate numerical solution of problems that arise from applications. The number of grid points used to discretize the computational domain has a significant impact on the speed of the optimization process. To take this issue into account, only 20 non-uniform cells are considered in our comparisons to test accuracy. The numerical results are compared in terms of integral moments, the average size of particles formed in the system and the number of particles in each cell. The integration of the discrete form of MFPT (2.3) is done using the ’MATLAB’ ’ODE45’ solver.

The average particle size (\(\bar{x}\)) can be calculated using the zeroth and first moments from the following relation:

$$\begin{aligned} \bar{x} = \displaystyle \frac{M_1(t)}{M_0(t)}. \end{aligned}$$
(4.1)

The absolute maximum relative error in the moments is calculated using the following expression:

$$\begin{aligned} \eta _i(t)=\left| \frac{M_i - \hat{M}_i}{M_i}\right| , \end{aligned}$$
(4.2)

where \(M_i\) and \(\hat{M}_i\) denote the exact and numerical ith moments, respectively. In addition, the sectional relative errors in the number density functions are estimated by means of the following relation:

$$\begin{aligned} \sigma _i(t)=\frac{\sum _{j=1}^{I} |N_j - \hat{N}_j| x_j^i}{\sum _{j=1}^{L} N_j x_j^i}. \end{aligned}$$
(4.3)

Moreover \(\sigma _i(t)\) provides the relative weighted sectional error in the number density function over the whole volume domain for \(i=0\).

Fig. 1
figure 1

Comparisons of numerical and exact results for binary breakage kernel with linear selection function using 20 cells

Full size image

4.1 Binary breakage kernel with linear selection function corresponding to exponential initial condition

For comparing the results, the binary breakage kernel with linear selection function (\(\Omega _t\left( x\right) =x\)) is chosen corresponding to an exponential initial condition \(n(0,x)=e^{-x}\). For running the numerical simulations, the computational size domain \(x_{\min }=10^{-9}\) to \(x_{\max }=3\) considered is subdivided into 20 non-uniform cells and the simulations are run till time \(t=5s\). One can see from Fig. 1a that the zeroth moment approximated by MFPT is more accurate than the cell average technique. Moreover, both numerical methods satisfy the mass conservation law, that is, the total mass remains consistent over the time domain (refer to Fig. 1b). In addition, the second order moment used to calculate the total area of the particles plays a significant role in many real life applications that arise in chemical engineering and pharmaceutical sciences is estimated with higher precision by the MFPT than the cell average technique (see Fig. 1c). The number of particles in each cell (\(N_i\)) calculated are plotted against the indices of the cell \(x_i\) in Fig. 1d. It shows that the MFPT and CAT exhibit almost equal accuracy but still the MFPT has the ability to trace the particles of larger size whereas the CAT fails to do so. However, it is worth noting that the the accuracy of the results can be improved to the desired level by assuming a more refined grid. The results obtained using a non-uniform grid with 40 cells are shown in Fig. 2. Still, the MFPT performs better than the existing method.

Fig. 2
figure 2

Comparisons of numerical and exact results for binary breakage kernel with linear selection function using 40 cells

Full size image
Table 1 Maximum error in different moments for binary breakage kernel with linear selection function for CAT and MFPT
Full size table

Moreover, the comparison is enhanced by calculating the relative errors (4.2) in the integral moments at the end time. From Table 1, it can be seen that the MFPT calculate these errors with more precision corresponding to fewer (20) non-uniform cells than the cell average technique. However, the two methods exhibit equal accuracy once the refined grid of 40 non-uniform cells is considered. In addition, to capture the errors in each cell of the computational domain, the sectional errors (4.3) are also estimated and listed in Table 2. Once again, the MFPT provides the numerical results with higher accuracy than the cell average technique.

4.2 Binary breakage kernel with quadratic selection function corresponding to monodisperse initial condition

In contrast to the previous cases, the binary breakage kernel with quadratic selection function (\(\Omega _t\left( x\right) =x^2\)) is chosen corresponding to the monodisperse initial condition \(n(0,x)=\delta (x-1)\). The computational grid and time required for running the simulation are taken to be same as in the previous case.

Table 2 Weighted error of number distribution for a binary breakage kernel with a linear selection function for the CAT and the MFPT
Full size table

One can observe from Fig. 3 that the MFPT performs better than the cell average technique in terms of calculating the zeroth order moment whereas other numerical results agree for the two methods. Similarly, the weighted sectional errors in the number density functions are calculated and listed in Table 3. Once again the MFPT shows better accuracy than the cell average technique for both coarse and refined grids consisting of 20 and 40 non-uniform grids, respectively.

It is worth mentioning here that for the other combination of binary breakage kernel and quadratic selection function corresponding to both exponential and monodisperse initial conditions, the numerical trends show similar behaviour for the case above. So the numerical results for these kernels are not presented here. However, to confirm the theoretical observations pertaining to the order of convergence, the experimental order of convergence is calculated for the binary breakage kernel with linear selection function corresponding to a monodisperse initial condition in the next section of the paper.

Fig. 3
figure 3

Comparisons of numerical and exact results for binary breakage kernel with quadratic selection function using 20 cells

Full size image

4.3 Experimental order of convergence

The theoretical results concerning the order of convergence of MFPT are tested against the numerical estimated order of convergence (EOC) for analytically tractable kernels similarly as by Kumar and Warnecke [8, 9] using the following formula:

$$\begin{aligned} {\text {EOC}} = \ln \left( \frac{E_I}{E_{2I}}\right) \Big /\ln (2), \end{aligned}$$
(4.4)

where \(E_I\) is the \(L^1\) error norm calculated by

$$\begin{aligned} E_I := \sum _{j=1}^{I}|N_j - \hat{N}_j |. \end{aligned}$$
(4.5)

The subscript I corresponds to the degrees of freedom, and the relative \(L^1\) error is computed by \(\displaystyle \left\| \textbf{N}- \hat{\textbf{N}} \right\| \Big / \Vert \textbf{N}\Vert \).

Table 3 Weighted error of number distribution for binary breakage kernel with quadratic selection function for the CAT and the MFPT
Full size table

The EOC is calculated for the binary breakage kernel with linear selection function corresponding to a monodisperse initial condition \(n(0,x)=\delta (x-1)\) for uniform and non-uniform grids (see Table 4). For calculating the EOC, the computational size domain \(x_{\min }=10^{-9}\) to \(x_{\max }=2\) considered is further subdivided into 30, 60, 120 and 240 non-uniform cells. All simulations are run until time \(t=5s\). Table 4 shows that the MFPT exhibits second order convergence on both uniform and non-uniforms grids similar to those consider by Kumar and Warnecke [8, 9].

Table 4 Experimental order of convergence corresponding to a monodisperse initial condition
Full size table

5 Concluding remarks and future prospects

In this work, the rate of convergence of the modified fixed pivot technique [11] is studied in detail. It is proved that the method has second order convergence rate irrespective of the spatial meshes used. Moreover, it has been also demonstrated that the rate of convergence is not dependent on the kind of breakage kernel and selection functions considered. Our numerical case study ascertains that the method exhibits improved accuracy compared to the existing cell average technique. This will allow practitioners to adapt the current approach for solving real life problems in the area of depolymerization, bubble column and granulation. Thanks to the mathematical flexibility of the modified fixed pivot technique in terms of not requiring specification of the total breakup rate of a mother particle and a daughter size distribution function at the start of the simulations, makes a strong case for the proposed method being easily coupled with CFD software.