1 Introduction

Fragmentation is a well-known mechanism occur in many applications of aerosol [10], depolymerization [1, 12], crystallization [31], pharmaceutical sciences [14, 15] and chemical engineering [23, 30]. The fragmentation process leads the formation of two or more smaller size particles after the disintegration of larger size particles. During this mechanism, particles of various sizes are formed in the system that can be tracked by the temporal change in the number density function governed by a mathematical population balance model. The total number of particles in the system increases during the fragmentation process whereas the total mass remains conserved.

The Mathematical model required to track the evolution of the number density function f(y,t) with time t > 0 corresponding to a fragmentation process is described by a linear integro-partial differential equation:

$$ \begin{array}{@{}rcl@{}} \frac{\partial f(y,t)}{\partial t}={\int}_{y}^{\infty} b(y,\xi) S(\xi) f(\xi,t) d\xi-S(y) f(y,t), \end{array} $$
(1)

subject to a initial condition

$$ \begin{array}{@{}rcl@{}} f(y,0)=f_{0}(y)\geq 0. \end{array} $$
(2)

Here, the selection function S(y) gives the fragmentation rate of the particles to be selected to undergo fragmentation mechanism. Moreover, the breakage kernel b(y,ξ) well known as fragmentation kernel describes the probability density function for the formation of particles having properties y from particles properties ξ. The first integral term on the right-hand side of (2) accounts for the birth of particle of size y due to the fragmentation of particle properties ξ and the second term represents the death of the particle size y due to fragmentation of that particle.

The breakage kernel b(y,ξ) must satisfies the following two properties:

$$ \begin{array}{@{}rcl@{}} {\int}_{0}^{\theta} b(y,\xi) dy&=v(\xi),~~ \text{for all}~y>0,~~~ b(y,\xi)=0,~~\text{for all}~y\geq\xi, \end{array} $$
(3)

and

$$ \begin{array}{@{}rcl@{}} {\int}_{0}^{\xi} yb(y,\xi) dy&=\xi,~ \forall~ \xi>0. \end{array} $$
(4)

The (3) is used to express the total number of fragments produced from a particles having size ξ and in general, v(ξ) ≥ 2. However, the relation (4) describes the mass conservation property, that is, when the particle of size ξ splits into smaller fragments then the total mass of those fragments is equal to ξ.

The main interest of this work is to approximate the number density function accurately. But for many real-life applications, specifically in chemical engineering and pharmaceutical sciences, the integral properties are of equal interest [33, 34, 37]. These integral properties are well known as moments calculated from the number density function by the following relation:

$$ \begin{array}{@{}rcl@{}} \mu_{i}(t) = {\int}_{0}^{\infty} y^{i} f(y,t) dy, \end{array} $$
(5)

where μi(t) gives the total number of particles (zeroth-order moment) for i = 0, however, for i = 1, the total mass in the system (first-order moment) can be calculated.

In this work, our main aim is to develop a new numerical approximation for the fragmentation (2). The other studies related to the existence [7, 25] and uniqueness [2, 24, 32] of the fragmentation equation are discussed in detail in these references. Despite of complex behavior of the fragmentation equation, some analytical solutions of fragmentation equation are derived by [16, 45, 46]. Other investigations concern scattering, self similarity and shattering are examined and discussed by [4, 5, 11]. In the literature, various researchers have developed different numerical methods to approximate the fragmentation (2). Among these numerical methods, the sectional methods [17, 22] are known for their accuracy to predict different order moments and the number density function. However, the complex mathematical formulations of these methods is their major drawback that makes these methods computationally very expensive.

Another simple class of approximations are method of moments [3, 27, 29, 44] in which the original equation is converted into a moment form and focused only on capturing the various order moments accurately. After the conversion of the original equation, the information related to the number density function is lost and other tools are required to track the number density function. In addition, stochastic methods [26, 43] are also very well known for their accuracy. But to achieve the accuracy, a large number of particles are required which make these methods computationally inefficient. However, the aforementioned issues were also addressed by developing highly efficient and accurate finite volume schemes [6, 8, 21, 38, 39, 41] which are highly efficient and accurate. Furthermore, these schemes are straightforward to extend to solve problems involving higher-dimensional population balances [9, 35, 36, 40, 42].

Recently, Kumar et al. [18] proposed a numerical scheme for solving a fragmentation equation by introducing weights into the discrete equation for achieving the mass conservation property and number preservation property (total number of particles in the system). However, the scheme does not capture the second-order moment accurately which plays very important role in estimating the total area of particles in many real-life applications [28]. Therefore, in the present work, we formulate a proposed finite volume scheme with a simpler mathematical formulation that is both computationally efficient and flexible to implement on uniform and nonuniform grids. The new method precisely approximates the zeroth- and first-order moments, as well as the number density function. Furthermore, the proposed scheme captures this particular moment with high accuracy without taking any measures for the accuracy of second-order moments. In contrast to the existing scheme [18], the proposed scheme is developed based on the notion of overlapping of cell [9].

Rest of the paper is structured as follows: in Section 2, the formulation of the proposed scheme is derived and the theoretical proofs of the mass conservation law and preservation of total number of particles are provided. In next Section 3, the convergence analysis of the proposed scheme is conducted on uniform and nonuniform grids. Section 4 is used to compare the accuracy of the proposed method against the analytical solutions for bench marking cases. In last Section 5, some important conclusions are made related to this study.

2 Numerical scheme

For developing the numerical method, the upper limit present in the first integral (2) must be replaced by a positive number say, ymax. Thus, the equation takes the following form:

$$ \begin{array}{@{}rcl@{}} \frac{\partial f(y,t)}{\partial t}={\int}_{y}^{y_{max}} b(y,\xi) S(\xi) f(\xi,t) d\xi-S(y) f(y,t). \end{array} $$
(6)

Further, it is necessary to discretize the continuous computational domain from 0 to ymax into L number of cells as shown in Fig. 1 where yL+ 1/2ymax. The lower and upper boundaries of the j th cell is denoted by yj− 1/2 and yj+ 1/2, respectively, for jN whose mean is \(y_{j} =\displaystyle \frac {y_{j-1/2}+y_{j+1/2}}{2}\). The step size of j th cell is given by Δyj = yj+ 1/2yj− 1/2.

Fig. 1
figure 1

One-dimensional domain discretization

Now let us discretize the time domain as tn+ 1 = tn + Δtn for all nN. Further assume that the average value of f at any time tn in j th cell is the approximation of the function \(f\left (y_{j},t^{n}\right )\). Further, the function f can be discretized as \(f_{j}(t)=f\left (y_{j},t\right )+\mathcal {O}\left ({\Delta } {y_{j}^{2}}\right )\) with the presumption that the function is sufficiently smooth. It is extremely difficult to solve this continuous (6) analytically due to the complexity of the original equation. As a result, the idea is to integrate the (6) over each domain of the cell j to transform the continuous equation into a set of ordinary differential equations, which further simplifies to

$$ \begin{array}{@{}rcl@{}} \frac{df_{j}}{dt}= {B_{j}}-{D_{j}},~ {\text {for}}~ j=1,2,\cdots,L. \end{array} $$
(7)

subject to a new initial condition

$$ \begin{array}{@{}rcl@{}} f_{j}(0)=\frac{1}{\Delta y_{j}} {\int}_{y_{j-1/2}}^{y_{j+1/2}} f(y,0)dy. \end{array} $$

Here the birth and death terms are written as

$$ \begin{array}{@{}rcl@{}} B_{j}=\frac{1}{\Delta y_{j}}{\int}_{y_{j-1/2}}^{y_{j+1/2}}{\int}_{y}^{y_{L+1/2}} b(y,\xi) S(\xi) f(\xi,t) d\xi dy, \end{array} $$
(8)

and

$$ \begin{array}{@{}rcl@{}} D_{j}=\frac{1}{\Delta y_{j}}{\int}_{y_{j-1/2}}^{y_{j+1/2}}S(y) f(y,t) dy. \end{array} $$
(9)

Consider the birth term (8) and simplify the equation by changing the order of integration, as shown below:

$$ \begin{array}{@{}rcl@{}} B_{j}=\frac{1}{\Delta y_{j}} \left[{\int}_{y_{j-1/2}}^{y_{j+1/2}} S(\xi) f(\xi){\int}_{y_{j-1/2}}^{\xi} b(y,\xi) dy d\xi +\sum\limits_{k=j+1}^{L} {\int}_{y_{k-1/2}}^{y_{k+1/2}} S(\xi) f(\xi){\int}_{y_{j-1/2}}^{y_{j+1/2}} b(y,\xi) dy d\xi \right]. \end{array} $$
(10)

The midpoint quadrature approximation is then employed to both integrals of (10) with respect to y, yielding:

$$ \begin{array}{@{}rcl@{}} B_{j}=S_{j} f_{j} {\int}_{y_{j-1/2}}^{y_{j}} b(y,y_{j})dy+\frac{1}{\Delta y_{j}}\sum\limits_{k=j+1}^{L}S_{k} f_{k}{\Delta} y_{k}{\int}_{y_{k-1/2}}^{y_{k+1/2}} b\left( y,y_{k}\right)dy +\mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$

In a similar manner, the discretize form of the death term (9) can be simplified to

$$ \begin{array}{@{}rcl@{}} D_{j}=S_{j} f_{j}+\mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$
(11)

Assume that in a j th cell, \(\hat {f}_{j}\) is the numerical approximation of the number density function fj. As a result, the discrete equations for approximating the number density function are as follows:

$$ \begin{array}{@{}rcl@{}} \frac{d\hat{f}_{j}}{dt} = \frac{1}{\Delta y_{j}}\sum\limits_{k=j}^{L} S_{k} \hat{f}_{k} {\Delta} y_{k} B_{jk} - S_{j} \hat{f}_{j}. \end{array} $$
(12)

Here

$$ \begin{array}{@{}rcl@{}} B_{jk}={\int}_{y_{j-1/2}}^{{p_{k}^{j}}} b\left( y,y_{k}\right) dy. \end{array} $$
(13)

where

$${p_{k}^{j}} = \begin{cases} y_{k}, & \text{when}~ k=j, \\ y_{k+1/2}, & \text{otherwise}. \end{cases}$$

Now we will define a proportionality constant in detail to introduce the notion of overlapping: For a uniform grid, the value of \({{\Lambda }_{j}^{k}}\) is 1 due to the fact that the newly formed particles fall completely inside any cell after the fragmentation mechanism. In contrast to a uniform grid, \({{\Lambda }_{j}^{k}}\) can be calculated based on the fact that when a particle having properties, say yk of a k th cell breaks, it form two particles having properties ykj and yj for a nonuniform grid. Then the lower and upper bounds of the newly formed particle after fragmentation mechanism becomes yj− 1/2 = yk− 1/2y(kj)− 1/2 and yj+ 1/2 = yk+ 1/2y(kj)+ 1/2. Practically, the possibility of formation of the particles completely inside the cell is rare. Hence, there are two possibilities of overlapping of the newly formed particle with one or more cells as given below:

  • When the upper boundary of the domain of the newly born particle will fall inside the cell and lower boundary is outside the cell, that is, yj− 1/2 > yk− 1/2y(kj)− 1/2 and yk+ 1/2y(kj)+ 1/2 < yj+ 1/2.

  • When the domain of newly born particle of size ykykj will totally fall inside the cell, that is, yj− 1/2yk− 1/2y(kj)− 1/2 and yk+ 1/2y(kj)+ 1/2yj+ 1/2.

The schematic representation of all possible cases of overlapping of the cells is shown in Fig. 2.

Fig. 2
figure 2

Representation of two basic possibilities for overlap

Mathematically, the proportionality constant for overlapping of the cells is calculated using the relation given below:

$$ \begin{array}{@{}rcl@{}} {{\Lambda}_{j}^{k}}=\frac{\overline{{\wedge^{k}_{j}}}-\underline{{\vee^{k}_{j}}}}{\Delta y_{k}-{\Delta} y_{k-j}}, \end{array} $$
(14)

where \(\overline {{\wedge ^{k}_{j}}}=\min \limits (y_{j+1/2}, y_{k+1/2}-y_{(k-j)+1/2})\) and \(\underline {{\vee ^{k}_{j}}}=\max \limits (y_{j-1/2},y_{k-1/2}-y_{(k-j)-1/2})\). Here, \(\overline {{\wedge ^{k}_{j}}}\) and \(\underline {{\vee ^{k}_{j}}}\) define the bounds of the intersection of the cells k and (kj) with cell j. The proportionality constant \({{\Lambda }_{j}^{k}}\) describes the extent of overlapping of the newly formed particle with cell j. The equality holds when the intersection is empty \(\left ({{\Lambda }_{j}^{k}}=0\right )\) and when the newly born particles domain falls completely inside the cell j \(\left ({{\Lambda }_{j}^{k}}=1\right )\).

Using the aforementioned relations and expressions, the formulation (12) becomes

$$ \begin{array}{@{}rcl@{}} \frac{d\hat{f}_{j}}{dt} = \frac{1}{\Delta y_{j}}\sum\limits_{k=j}^{L} S_{k} \hat{f}_{k} {\Delta} y_{k} B_{jk}{{\Lambda}_{j}^{k}} - S_{j} \hat{f}_{j}. \end{array} $$
(15)

It can be noticed that the (15) does not give any account for the mass conservation law which is a necessary condition for any numerical method. This signifies that the above formulation cannot be used to track the true behavior of the number density and require some adjustments to predict the integral properties (zeroth- and first-order moments) accurately corresponding to the number density function. This issue will be resolved by adding two weights in the birth and death terms of the (15) which gives

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \hat{f}_{j}= \underbrace{\frac{1}{\Delta y_{j}}\sum\limits_{k=j}^{L} S_{k} \hat{f}_{k} B_{jk} {{\Lambda}_{j}^{k}} {\Delta} y_{k} {\omega_{k}^{b}}}_{{\widehat{B}_{j}}}- \underbrace{\frac{1}{\Delta y_{j}}S_{j} \hat{f}_{j} {\Delta} y_{j} {\omega_{j}^{d}}}_{{\widehat{D}_{j}}} . \end{array} $$
(16)

Here the weights responsible for the preservation of total number of particles and conservation of the total mass in the system are given by

$$ \begin{array}{@{}rcl@{}} {\omega_{k}^{b}}=\frac{1}{{{\Lambda}_{j}^{k}}}\frac{y_{k}(v(y_{k})-1)}{{\sum}_{j=1}^{k-1}\left( y_{k}-y_{j}\right)B_{jk}}, \end{array} $$
(17)

and

$$ \begin{array}{@{}rcl@{}} {\omega_{k}^{d}}=\frac{{\omega_{k}^{b}}}{y_{k}}\sum\limits_{j=1}^{k}y_{j}B_{jk}, ~~j=2,3,\cdots, L. \end{array} $$
(18)

The values of weights \({\omega _{1}^{b}}\) and \({\omega _{1}^{d}}\) are consider to be zero.

Using the above notations in the (16), the final expression for a proposed finite volume scheme for solving fragmentation equation is given by

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \hat{f}_{j}= \underbrace{\frac{1}{\Delta y_{j}}\sum\limits_{k=j}^{L} S_{k} \hat{f}_{k} B_{jk} {{\Lambda}_{j}^{k}} {\Delta} y_{k} {\omega_{k}^{b}}}_{{\widehat{B}_{j}}}- \underbrace{\frac{1}{\Delta y_{j}}S_{j} \hat{f}_{j} {\Delta} y_{j} {\omega_{k}^{d}}}_{{\widehat{D}_{j}}} . \end{array} $$
(19)

The theoretical proofs of the mass conservation property and number preservation property are given below. The numerical scheme holds the number preservation property if it satisfies the following relation:

$$ \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} {\Delta} y_{j}=\sum\limits_{j=1}^{L}S_{j} \hat{f}_{j}{\Delta} y_{j}(v(y_{j})-1). $$
(20)

However, the numerical scheme holds the mass conservation property when it satisfies the following condition:

$$ \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} y_{j}{\Delta} y_{j}=0. $$
(21)

Equations (20) and (21) can be easily obtained from (7) by multiplying with yj and taking summation corresponding to j = 0 and j = 1, respectively.

Proposition 1

The discrete formulation (19) holds the number preservation property (zeroth-order moment).

Proof

Take summation over all j on the discrete formulation (19), we get

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} {\Delta} y_{j} =\sum\limits_{j=1}^{k} \sum\limits_{k=j}^{L} S_{k} \hat{f}_{k} {\Delta} y_{k} {{{\Lambda}_{j}^{k}}} B_{jk}{\omega_{k}^{b}} -\sum\limits_{j=1}^{I} S_{j} \hat{f}_{j} {\Delta} y_{j}{\omega_{k}^{d}}. \end{array} $$
(22)

Now changing the order of the summation, we have

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} {\Delta} y_{j} =\sum\limits_{k=1}^{L} \sum\limits_{j=1}^{k} S_{k} \hat{f}_{k} {\Delta} y_{k} {{{\Lambda}_{j}^{k}}} B_{jk}{\omega_{k}^{b}} -\sum\limits_{j=1}^{I} S_{j} \hat{f}_{j} {\Delta} y_{j}{\omega_{k}^{d}}. \end{array} $$
(23)

Substituting the values of weight from (18), the above equation gives

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} {\Delta} y_{j}= \sum\limits_{k=1}^{L} S_{k} \hat{f}_{k} {\Delta} y_{k}\left[ \sum\limits_{j=1}^{k} B_{jk}{\omega_{k}^{b}}- \frac{{\omega_{k}^{b}}}{y_{k}}\sum\limits_{j=1}^{k}y_{j}B_{jk}\right]. \end{array} $$

Further replace the values of the weight from (17) and after simplification, it gives

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} {\Delta} y_{j}= \sum\limits_{k=1}^{L} S_{k} \hat{f}_{k} {{\Lambda}_{j}^{k}}{\Delta} y_{k} \frac{y_{k}[v(y_{k})-1]}{y_{k}{{\Lambda}_{j}^{k}}{\sum}_{k=1}^{j}\left( y_{k}-y_{j}\right)B_{jk}}\left( \sum\limits_{j=1}^{k}\left( y_{k}-y_{j}\right)B_{jk}\right). \end{array} $$

This implies

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} \hat{f}_{j} {\Delta} y_{j}=\sum\limits_{k=1}^{L}S_{k} \hat{f}_{k}{\Delta} y_{k}(v(y_{k})-1). \end{array} $$
(24)

This shows that the proposed finite volume scheme is preserving the zeroth-order moment. □

Proposition 2

The finite volume scheme (19) is conserving the total mass in the system, that is, first-order moment.

Proof

For proving the mass conservation property, the (19) should be multiplied with yj on both sides and take summation over j will give

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} y_{j}\hat{f}_{j} {\Delta} y_{j} =\sum\limits_{k=1}^{I} \sum\limits_{j=1}^{k} S_{k} \hat{f}_{k} {\Delta} y_{k} {{{\Lambda}_{j}^{k}}} B_{jk}{\omega_{k}^{b}} -\sum\limits_{j=1}^{I} S_{j} \hat{f}_{j} {\Delta} y_{j}{\omega_{k}^{d}}. \end{array} $$
(25)

The above can also be rewritten as

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} y_{j}\hat{f}_{j} {\Delta} y_{j} =\sum\limits_{k=1}^{I} S_{k} \hat{f}_{k} {\Delta} y_{k} {{{\Lambda}_{j}^{k}}} \left( {\omega_{k}^{b}}\sum\limits_{j=1}^{k}B_{jk} -{\omega_{k}^{d}} x_{k}\right). \end{array} $$
(26)

Substituting the values of \({\omega _{k}^{d}}\) from (18) and after simplification, we get

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{j=1}^{L} y_{j} \hat{f}_{j} {\Delta} y_{j} =\sum\limits_{k=1}^{I} S_{k} \hat{f}_{k} {\Delta} y_{k} {{{\Lambda}_{j}^{k}}}{\omega_{k}^{b}}\left( \sum\limits_{j=1}^{k} x_{j}B_{jk} -\sum\limits_{j=1}^{k}x_{j}B_{jk}\right). \end{array} $$
(27)

This further gives

$$ \begin{array}{@{}rcl@{}} \frac{d}{dt} \sum\limits_{i=1}^{L} y_{j}\hat{f}_{j}{\Delta} y_{j}=0. \end{array} $$
(28)

This implies that the total mass in the system is conserved, that is, no mass is lost from the system. □

In the next section of the article, the investigation of convergence analysis of the numerical scheme will be carried out.

3 Convergence analysis

It is necessary to write the equations in vector form before discussing the convergence analysis of the proposed finite volume scheme. Assume that the vectors f and \(\hat {\mathbf {f}}\) denote, respectively, the average exact and average numerical values of the number density functions. In vector form, the discrete (19) can be written as

$$ \begin{array}{@{}rcl@{}} \frac{\partial \hat{\mathbf{f}}}{\partial t} = \mathbf{J}(\hat{\mathbf{f}}),~~~\hat{\mathbf{f}}(0)= \mathbf{f}(0). \end{array} $$
(29)

Here, the birth term \((\hat {\mathbf {B}})\), death term \((\hat {\mathbf {D}})\) and \(\mathbf {J} \in \mathcal {R}^{I}\) are the functions of \(\hat {f}\) with the components

$$ \begin{array}{@{}rcl@{}} \hat{B}_{j}(\hat{\mathbf{f}})= \frac{1}{ {\Delta} y_{j}} \sum\limits_{k=i}^{L} {\omega_{k}^{b}} S_{k} \hat{f}_{k} {\Delta} y_{k} {B}_{jk}, \end{array} $$
(30)

and

$$ \begin{array}{@{}rcl@{}} \hat{D}_{j}(\hat{\mathbf{f}})= {\omega_{j}^{d}} S_{j} \hat{f}_{j}. \end{array} $$
(31)

Therefore, the final form of equation is

$$ \begin{array}{@{}rcl@{}} {J_{j}}(\hat{\mathbf{f}})= \hat{B}_{j}\left( \hat{f}\right)-\hat{D}_{j}\left( \hat{f}\right). \end{array} $$
(32)

To conduct the convergence of the discrete system, first it is necessary to define the norm L1 which is given by

$$ \begin{array}{@{}rcl@{}} \Vert \mathbf{f}(t)\Vert = \sum\limits_{j=1}^{L}|f_{j}(t)|{\Delta} y_{j}. \end{array} $$
(33)

Now, let us give some important definitions which will be required in conducting the convergence analysis of the proposed finite volume scheme.

Definition 1

Spatial Truncation Error σ(t) is obtained by substituting the exact solution \(\textbf {f}=[f_{1}(t),f_{2}(t),\dots ,f_{L}(t)]\) in the discrete system of equations, that is,

$$ \begin{array}{@{}rcl@{}} \sigma(t)=\frac{d\mathbf{f}(t)}{dt}- \mathbf{J}(\mathbf{f}). \end{array} $$

The numerical scheme is said to be consistent of order p, if \({\Delta } y \rightarrow 0\)

$$ \begin{array}{@{}rcl@{}} \Vert \sigma(t) \Vert=\mathcal{O} ({\Delta} y^{p}),~~~\text{uniformly for all } t,~ 0\leq t\leq T. \end{array} $$

Now, let us define another type of discretization error which will be used to find the order of convergence.

Definition 2

Global Discretization Error for any numerical scheme is the difference between the exact and numerical solution \(\epsilon (t)=f(t)-\hat {f}(t)\). The numerical scheme is said to be convergent of order p if, for \({\Delta } y\rightarrow 0\),

$$ \begin{array}{@{}rcl@{}} \Vert \epsilon(t) \Vert=\mathcal{O} ({\Delta} y^{p}),~~~\text{uniformly for all } t,~ 0\leq t\leq T. \end{array} $$
(34)

Proposition 3

Let us suppose that a Lipschitz condition on J(f) is satisfied for 0 ≤ tT and for all \(f,\hat {f}\in \mathbb {R}^{L}\). That is, J satisfies

$$ \begin{array}{@{}rcl@{}} \Vert\boldsymbol{J(f)}-\boldsymbol{J}(\hat{\boldsymbol{f}})\Vert \leq M\Vert \boldsymbol{f}-\hat{\boldsymbol{f}} \Vert, ~~{M} < \infty. \end{array} $$

Then a consistent discretization scheme is also convergent and the convergent order is same as the order of consistency.

First, it is important to prove the Lipschitz condition for the function J.

Proposition 4

Let us assume that the kernels S and b are twice continuously differentiable functions over ]0,ymax] and ]0,ymax], respectively, then there exist a constant

$$ M=C \underset{y \in (0,y_{max}]}{\max} S(y) v(y)< \infty,~~\text{where}~~ C>0, $$

such that the Lipschitz condition on J is satisfied for all \(\boldsymbol {f}~ and ~\hat {\boldsymbol {f}} \in \mathbb {R}^{L}\), that is,

$$ \begin{array}{@{}rcl@{}} \Vert \boldsymbol{J(f)}-\boldsymbol{J} (\hat{\boldsymbol{f}}) \Vert \leq M \Vert {\boldsymbol{f}}-\hat{\boldsymbol{f}} \Vert. \end{array} $$

Proof

Further consider the definition of the norm defined in (33)

$$ \begin{array}{@{}rcl@{}} \Vert\mathbf{J(f)}-\mathbf{J}(\hat{\mathbf{f}})\Vert =\sum\limits_{j=1}^{L}|J_{j}(f_{j})- J_{j}\left( \hat{f}_{j}\right)| \end{array} $$

This implies

$$ \Vert\mathbf{J(f)}-\mathbf{J}(\hat{\mathbf{f}})\Vert =\sum\limits_{j=1}^{L} \frac{1}{\Delta y_{j}}\left|\sum\limits_{k=j}^{L} S_{k} B_{jk} {\Delta} y_{k} {{\Lambda}_{j}^{k}} \left( f_{k}-\hat{f}_{k}\right){\omega_{k}^{b}}-S_{j} \left( f_{j}-\hat{f}_{j}\right) {\Delta} y_{j} {\omega_{j}^{d}}\right|{\Delta} y_{j}, $$

which further simplifies to

$$ \begin{array}{@{}rcl@{}} \Vert\mathbf{J(f)}-\mathbf{J}\left( \hat{\mathbf{g}}\right)\Vert \leq \underbrace {\sum\limits_{j=1}^{L} \sum\limits_{k=j}^{L} S_{k} B_{jk} {\Delta} y_{k} {{\Lambda}_{j}^{k}} |f_{k}-\hat{f}_{k}| {\omega_{k}^{b}}}_{\text{$T_{1}$}} + \underbrace{\sum\limits_{j=1}^{L} S_{j} |f_{j}-\hat{f}_{j}| {\Delta} y_{j} {\omega_{j}^{d}}}_{\text{$T_{2}$}}. \end{array} $$
(35)

Now simplify the first term by changing the order of the summation, we get

$$T_{1} = \sum\limits_{j=1}^{L} \sum\limits_{k=j}^{L} S_{k} B_{jk} {\Delta} y_{k} {{\Lambda}_{j}^{k}} |f_{k}-\hat{f}_{k}| {\omega_{j}^{b}}= \sum\limits_{k=1}^{L} \sum\limits_{j=1}^{k} S_{k} B_{jk} {\Delta} y_{k} {{\Lambda}_{j}^{k}} |f_{k}-\hat{f}_{k}| {\omega_{j}^{b}}. $$

After substituting the value of \({\omega _{k}^{b}}\) and using the fact that \({\omega _{1}^{b}}=0\) in the above equation, we have

$$ \begin{array}{@{}rcl@{}} T_{1} &=& \sum\limits_{k=2}^{L} \sum\limits_{k=i}^{L} S_{k} B_{ik} {{\Lambda}_{j}^{k}}|f_{k}-\hat{f}_{k}| {\Delta} y_{k} \frac{1}{{{\Lambda}_{j}^{k}}}\left[ \frac{y_{k}v(y_{k})}{ {\sum}_{j=1}^{k} \left( y_{k}-y_{j}\right)B_{jk}}-\frac{y_{k}}{ {\sum}_{j=1}^{k} \left( y_{k}-y_{j}\right)B_{jk}} \right], \\ &=& \sum\limits_{k=2}^{L} \sum\limits_{k=i}^{L} S_{k} B_{ik} |f_{k}-\hat{f}_{k}| {\Delta} y_{k}\left[ \frac{\left( y_{k}-y_{j}\right)v(y_{k})}{ {\sum}_{j=1}^{k} \left( y_{k}-y_{j}\right)B_{jk}}-\frac{y_{k} v(y_{k})}{ {\sum}_{j=1}^{k} \left( y_{k}-y_{j}\right)B_{jk}} \right] \\&& -\sum\limits_{k=2}^{L} \sum\limits_{k=i}^{L} S_{k} B_{ik}|f_{k}-\hat{f}_{k}| {\Delta} y_{k}\frac{y_{k}}{ {\sum}_{j=1}^{k} \left( y_{k}-y_{j}\right)B_{jk}}. \end{array} $$
(36)

Rearranging the above equation, we get

$$ \begin{array}{@{}rcl@{}} T_{1} \leq & \sum\limits_{k=2}^{L} v(y_{k})S_{k} |f_{k}-\hat{f}_{k}| {\Delta} y_{k} + \sum\limits_{k=2}^{L} \frac{v(y_{k})S_{k} B_{ik} |f_{k}-\hat{f}_{k}| {\Delta} y_{k}}{ {\sum}_{j=1}^{k} \left( y_{k}-y_{j}\right)B_{jk}} \left[ \sum\limits_{j=1}^{k} y_{j}B_{jk} -y_{k}\right], \end{array} $$
(37)

As the breakage function b(y,ξ) is considered to be twice continuously differentiable function and further using the midpoint and the right end approximation of the integrals, we obtain

$$ \begin{array}{@{}rcl@{}} y_{k}={\int}_{0}^{y_{k}} b(y, y_{k})dy = \sum\limits_{j=1}^{k} {\int}_{y_{j-1/2}}^{{p_{k}^{j}}} yb(y, y_{k})dy= \sum\limits_{j=1}^{k} y_{j} B_{jk}+\mathcal{O} ({\Delta} y)^{2}. \end{array} $$

Hence,

$$ \begin{array}{@{}rcl@{}} \left|y_{j} - \sum\limits_{j=1}^{k} y_{j} B_{jk}\right| = \mathcal{O} ({\Delta} y)^{2} \leq M_{1} ({\Delta} y)^{2}, ~\text{where} ~M_{1}<\infty ~\text{is a constant}. \end{array} $$
(38)

Additionally, for k = 2,3,⋯ ,L, we have

$$ \begin{array}{@{}rcl@{}} \sum\limits_{j=1}^{k} \left( y_{k}-y_{j}\right) B_{jk} \geq (y_{k}-y_{k-1}) \sum\limits_{j=1}^{k-1} B_{jk} \geq {\Delta} y_{k} {\int}_{0}^{y_{k-1/2}} b\left( y,y_{k}\right)dy \geq M_{2} (\delta x) \end{array} $$
(39)

where M2 is a constant, satisfying

$$ \begin{array}{@{}rcl@{}} 0<M_{2}= \underset{k\in \{2,3,\cdots, L\}}{\min} \left[ {\int}_{0}^{y_{k-1/2}} b\left( y,y_{k}\right)dy \right] \end{array} $$
(40)

Using the above relations in (37), we get

$$ \begin{array}{@{}rcl@{}} T_{1} &\leq& \sum\limits_{k=2}^{L} v(y_{k})S_{k} |f_{k}-\hat{f}_{k}| {\Delta} y_{k} + \sum\limits_{k=2}^{L} \frac {({\Delta} y)^{2}} {\delta y}\frac {M_{1}} {M_{2}} v(y_{k})S_{k} |f_{k}-\hat{f}_{k}| {\Delta} y_{k} \\& \leq& \left[ 1+ \frac {M_{1}} {M_{2}} K y_{max}\right] \sum\limits_{k=1}^{L} v(y_{k})S_{k} |f_{k}-\hat{f}_{k}| {\Delta} y_{k} \end{array} $$
(41)
$$ \begin{array}{@{}rcl@{}}& \leq &\left[ 1+ \frac {M_{1}} {M_{2}} K y_{max}\right] \underset{{k\in \{2,3,\cdots, L\}}}{\max} [v(y)S(y)] \sum\limits_{k=1}^{L}|f_{k}-\hat{f}_{k}| {\Delta} y_{k} \end{array} $$
(42)

Here \(\delta y =\min \limits _{i} {\Delta } y_{i}\). The above relation can further be simplified to

$$ \begin{array}{@{}rcl@{}} T_{1} \leq M_{3} ||\mathbf{f} -\hat{\mathbf{f}}||. \end{array} $$
(43)

where \(M_{3} = \frac {2}{W} \max \limits _{\substack {y\in \{2,3,\cdots , y_{max}\}}} [v(y)S(y)] \leq \infty \) and \(W = 2\left [1+\frac {M_{1}} {M_{2}} K y_{max}\right ]\).

Further, consider the second term T2 from the (35),

$$ \begin{array}{@{}rcl@{}} T_{2} = \sum\limits_{j=1}^{L} S_{j} |f_{j}-\hat{f}_{j}| {\Delta} y_{j} {\omega_{k}^{d}} = \sum\limits_{j=1}^{L} S_{j} |f_{j}-\hat{f}_{j}| {\Delta} y_{j} y_{k} \frac{{\omega_{k}^{b}}}{y_{k}}\sum\limits_{j=1}^{k}B_{jk}. \end{array} $$

Using the fact that yjykj = 1,2,3,⋯ ,k,

$$ \begin{array}{@{}rcl@{}} T_{2} \leq \sum\limits_{k=1}^{L} S_{j} |f_{j}-\hat{f}_{j}| {\Delta} y_{j} \frac{{\omega_{k}^{b}}}{y_{k}}\sum\limits_{j=1}^{k}y_{j}B_{jk}. \end{array} $$

Now changing the order of the integration leads to

$$ \begin{array}{@{}rcl@{}} T_{2} \leq \sum\limits_{j=1}^{L}\sum\limits_{k=j}^{L} S_{j} |f_{j}-\hat{f}_{j}| {\Delta} y_{j} \frac{{\omega_{k}^{b}}}{y_{k}}y_{j}B_{jk} =T_{1}. \end{array} $$

This gives

$$ \begin{array}{@{}rcl@{}} T_{2} \leq M_{3} ||\mathbf{f} -\hat{\mathbf{f}}||. \end{array} $$
(44)

Using the (43) and (44), the (35) takes the following form:

$$ \begin{array}{@{}rcl@{}} \Vert\mathbf{J(f)}-\mathbf{J}\left( \hat{\mathbf{g}}\right)\Vert \leq M ||\mathbf{f} -\hat {\mathbf{f}}||, \end{array} $$
(45)

where \(M=2M_{3} < \infty \) is a Lipschitz constant. □

Now, the next aim is to prove the order of convergence of the proposed finite volume scheme by stating the following theorem.

Theorem 1

Suppose that the functions S and b are twice continuously differentiable functions over (0,ymax] and (0,ymax] × (0,ymax], respectively. Then, the solution of the discretization (29) is non-negative and consistent, with a truncation error of order 2, independent of the type of grid. Moreover, the scheme is convergent and the order of convergence is the same as the order of consistency.

Proof

For proving the theorem, we need to prove the three components of the solution, that is, non-negativity, consistency, and convergence. First we will begin the exercise with the non-negativity of the solution. □

Non-negativity

For any non-negative number density function \(\hat {f} \in \mathbb {R}^{L}\), (for all \(\hat {f} \geq 0\) which has j th component equals to zero). So, (30) and (31) give

$$ \begin{array}{@{}rcl@{}} \hat{B}_{j}(\hat{\mathbf{f}}) \geq 0~~~\text{and}~~~\hat{D}_{j}(\hat{\mathbf{f}}) = 0. \end{array} $$

Then, the (32) implies \({J_{j}}\left (\hat { \mathbf {f}}\right )= \hat {B}_{j}\left (\hat {\mathbf {f}}\right )-\hat {D}_{j}(\hat {\mathbf {f}})\geq 0\). For all \(i=1, 2,\dots , L\), using the Theorem 7.1, Chapter 1 of [13] and Proposition 3 give the non-negativity of the solution.

Consistency

Using the Definition 1, the j th component of the spatial truncation error can be written as

$$ \begin{array}{@{}rcl@{}} \sigma_{j}(t) = \frac{df_{j}(t)}{dt}-J_{j}(f_{j}(t)). \end{array} $$

Further using (7) and (32), the above relation becomes

$$ \begin{array}{@{}rcl@{}} \sigma_{i}(t) = \underbrace{B_{j}-\widehat{B}_{j}}_{I_{1}}-\underbrace{(D_{j}-\widehat{D}_{j})}_{I_{2}}. \end{array} $$
(46)

Now consider the first term I1 of the (46),

$$ \begin{array}{@{}rcl@{}} B_{j}-\widehat{B}_{j} = \frac{1}{\Delta y_{j}} \left[\sum\limits_{k=j}^{L} S_{k} f_{k} {\Delta} y_{k} B_{jk} - \sum\limits_{k=j}^{L} S_{k} f_{k} {\Delta} y_{k} B_{jk} {{\Lambda}_{j}^{k}} {\omega_{k}^{b}} \right] + \mathcal{O}\left( {\Delta} y^{2}\right), \end{array} $$
(47)

Combining the terms, we get

$$ \begin{array}{@{}rcl@{}} B_{j}-\widehat{B}_{j} = \frac{1}{\Delta y_{j}} \left[\sum\limits_{k=j}^{L} S_{k} f_{k} {\Delta} y_{k} B_{jk} \left( 1-{{\Lambda}_{j}^{k}}{\omega_{k}^{b}}\right) \right] + \mathcal{O}\left( {\Delta} y^{2}\right), \end{array} $$
(48)

Further estimate the order of the term \(1-{{\Lambda }_{j}^{k}}{\omega _{k}^{b}}\) as follows:

$$ \begin{array}{@{}rcl@{}} 1-{{\Lambda}_{j}^{k}}{\omega_{k}^{b}} = {{\Lambda}_{j}^{k}} \frac{{y_{k}\left( {\sum}_{j=1}^{k}B_{jk} -1\right)}}{{{{\Lambda}_{j}^{k}}}{{\sum}_{j=1}^{k}\left( y_{k}-y_{j}\right)B_{jk}}} = \frac{y_{k} - {\sum}_{j=1}^{k}y_{j}B_{jk}}{{\sum}_{j=1}^{k}\left( y_{k}-y_{j}\right)B_{jk}}. \end{array} $$

Using the mass conservation property of the breakage function defined in (3) gives

$$ \begin{array}{@{}rcl@{}} 1-{{\Lambda}_{j}^{k}}{\omega_{k}^{b}} = \frac{{\sum}_{j=1}^{k}{\int}_{y_{1/2}}^{{p_{k}^{j}}} \left( y-y_{j}\right)b\left( y,y_{k}\right)dy}{{\sum}_{j=1}^{k}\left( y_{k}-y_{j}\right)B_{jk}} \end{array} $$
(49)

Apply the midpoint and right end quadrature approximations, it can be easily shown that the numerator of above equation is of order 2 whereas the denominator exhibits order 0. Hence, we have

$$ \begin{array}{@{}rcl@{}} 1-{{\Lambda}_{j}^{k}}{\omega_{k}^{b}} = \mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$

Using the above relation, the (48) implies

$$ \begin{array}{@{}rcl@{}} B_{j}-\widehat{B}_{j} = \mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$
(50)

Similar to the term I1, now let us discuss the order of consistency for the term I2 of (46)

$$ \begin{array}{@{}rcl@{}} D_{j}-\widehat{D}_{j}=\left( 1-{\omega_{j}^{d}}\right)S_{j} f_{j}+\mathcal{O}\left( {\Delta} y^{2}\right) \end{array} $$
(51)

Substituting the value of \({\omega _{j}^{d}}\) from (18) in the above equation, we get

$$ \begin{array}{@{}rcl@{}} D_{j}-\widehat{D}_{j}&=&\left( 1-\frac{{\omega_{j}^{d}}}{y_{j}}\sum\limits_{j=1}^{k}y_{j}B_{j,k}\right) S_{j} f_{j}+\mathcal{O}\left( {\Delta} y^{2}\right), \end{array} $$
(52)
$$ \begin{array}{@{}rcl@{}} &=&\frac{1}{y_{j}}\left( y_{j}-{\omega_{j}^{b}}\sum\limits_{j=1}^{k}y_{j}B_{jk}\right) S_{j} f_{j}+\mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$
(53)

It can be noted that

$$ \begin{array}{@{}rcl@{}} y_{j}={\int}_{0}^{y_{j}} y b(y,y_{j}) dm= \sum\limits_{j=1}^{k} {\int}_{y_{j-1/2}}^{{p_{j}^{k}}} y b\left( y,y_{k}\right) dy = \sum\limits_{k=1}^{j} y_{j} B_{kj}+\mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$

Using yj from the above equation in (51) will give

$$ \begin{array}{@{}rcl@{}} D_{j}-\widehat{D}_{j}=\left( 1-{\omega_{j}^{b}}\right)\left( \frac{1}{y_{j}}\sum\limits_{j=1}^{k}B_{jk}\right) S_{j} f_{j}+\mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$
(54)

As proved earlier that \(1-{\omega _{j}^{b}} = \mathcal {O}\left ({\Delta } y^{2}\right )\), we have

$$ \begin{array}{@{}rcl@{}} D_{j}-\widehat{D}_{j}=\mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$
(55)

Hence, from (46), (50) and (55), we get

$$ \begin{array}{@{}rcl@{}} \sigma_{j}(t)=\mathcal{O}\left( {\Delta} y^{2}\right) \Rightarrow \Vert {\sigma}(t) \Vert = \mathcal{O}\left( {\Delta} y^{2}\right). \end{array} $$

Convergence

From Propositions 3 and 4, and the above result on consistency proves that the proposed finite volume scheme converges to the order same as the order of consistency.

In the next section, the numerical results for various combinations of selection functions and breakage kernels will be discussed.

4 Numerical validation

To validate the proposed finite volume scheme, the numerical results in terms of different order moments as well as the number density functions are compared against several benchmark cases using nonuniform grids. In particular two different cases will be tested (a) analytically tractable kernels for which the analytical results for both moments and number density function are available in the literature, and (b) practically oriented kernel for which analytical results are not available. The analytical solutions of moment and number density functions corresponding to the different initial conditions are available in literature [16, 45, 46]. For our comparison, monodisperse f(y,0) = δ(y − 1) initial conditions are considered for analytically tractable cases, whereas for the practically oriented problem, the following initial condition is considered:

$$ \begin{array}{@{}rcl@{}} f(y,0) = \frac{1}{\sigma \sqrt2\pi} \exp\left[ \frac{-(y-\mu)^{2}}{2\sigma^{2}}\right], ~\sigma^{2} =0.01, ~ \mu =0.5. \end{array} $$
(55)

The comparison is enhanced by quantifying the weighted sectional errors exist in the number density function which can be estimated using the following relation:

$$ \begin{array}{@{}rcl@{}} \eta_{i}(t)=\frac{{\sum}_{j=1}^{L} |f_{j}^{exc}-f_{j}^{num}| {y_{j}^{i}} {\Delta} {y_{j}^{i}}}{{\sum}_{j=1}^{L} f_{j}^{exc} {y_{j}^{i}} {\Delta} {y_{j}^{i}} }, \end{array} $$
(55)

where σi(t) describes the relative weighted sectional error in the number density function over the whole volume domain for i = 0. Similarly, other order relative sectional errors in the number density function can be defined. These errors are evaluated for those cases whose analytical solutions are available in the literature and calculated at the end time. The integration of discrete form of breakage population balance (19) is solved using MATLAB ODE15s solver.

Moreover, the theoretical aspect of the convergence of the proposed finite volume scheme is also tested numerically by calculating the Experimental Order of Convergence (EOC) for analytically tractable kernels using the following expression given by [19, 20]:

$$ \begin{array}{@{}rcl@{}} {\text{EOC}} = {\ln}\left( \frac{E_{L}}{E_{2L}}\right)/{\ln}(2), \end{array} $$
(55)

where EI and E2I describe the L1 error norm calculated by

$$ \begin{array}{@{}rcl@{}} \sum\limits_{j=1}^{L}|f_{j}^{exc}-f_{j}^{num} |{\Delta} {y_{j}^{i}}. \end{array} $$
(55)

Here, \(f_{j}^{exc}\) and \( f_{j}^{num}\) describe the number of particles obtained analytically and numerically, respectively, and the symbols L and 2L correspond to the number of degrees of freedom.

4.1 Test case I

Let us begin the testing of the proposed finite volume scheme for the binary breakage kernel b(y,ξ) = 2/ξ and linear selection function S(y) = y corresponding to the monodisperse initial condition. The computational domain considered to run the numerical simulation begins with ymin = 10− 9 and ends with ymax = 1 is divided into 30 nonuniform cells. The advantages of the nonuniform grids over uniform grids is provided in detail by [9]. It has been shown that a large number of grid points are required for the uniform grid to capture the integral properties accurately. The simulations are run between time 0 s and 1000 s where the extent of breakage attained is 1047.91, that is, \(\displaystyle \frac {\mu (t)}{\mu (0)}\approx 1047.91\), with μ(t) denotes the zeroth-order moment at time t. It is important to note that the numerical scheme developed by [9] for solving aggregation population balance equation stopped their simulations when the time reached t = 1.5 s.

The qualitative comparison of numerical results in terms of number density function and different order moments against its analytical results is demonstrated in Fig. 3. It shows that the zeroth- and first-order moments are well predicted by the proposed scheme and matching well with the analytical results, that is, the proposed scheme preserves the total number of the particles as well as conserves the total mass in the system (see Fig. 3a). Moreover, the second-order moment computed by the proposed scheme also overlaps with the analytical result as shown in Fig. 3b. It is important to notice that the proposed method accounts the accuracy of the zeroth- and first-order moments. However, the proposed finite volume scheme still predicts the second-order moment very accurately. Additionally, the numerical scheme developed by [18] focused only on zeroth- and first-order moments, but did not give any account for the accuracy of the second-order moment. This shows that the proposed approximation has the ability to track the higher order moments accurately. Also, the number density function plotted in Fig. 3c reveals that the finite volume scheme is in excellent agreement with the analytical result.

Fig. 3
figure 3

Comparison of various order moments and number density for binary breakage kernel, i.e., b(m,m) = 2/m with linear selection function S(m) = m for monodisperse initial condition

In order to quantify the errors in the number density function, the weighted sectional errors calculated using the expression (4) are listed in Table 1 for two different size grids. It exhibits that these errors decrease to 50% when a refined grid with 60 nonuniform cells is used. To ensure the order of the convergence of the proposed scheme, the comparison of the numerical experimental order of convergence is provided for both uniform and nonuniform grids in Table 2. It illustrates that the proposed scheme shows second-order convergence irrespective of the kind of grids used for discretization. It can also be seen from the table that the error obtained using L1 norm (4) decreases enormously as the number of cells increased in the domain.

Table 1 Weighted error of number density for binary breakage kernel with quadratic selection function
Table 2 Experimental order of convergence for Case I corresponding to monodisperse initial condition

4.2 Test case II

Similar to the previous case, the comparison of the numerical moments and number density function with analytical results is conducted for binary breakage kernel and quadratic selection function S(y) = y2. The computational size domain taken from ymin = 10− 9 and ymax = 1.5 is partitioned into 30 nonuniform cells. The extent of breakage achieved is \(\displaystyle \frac {\mu (t)}{\mu (0)}\approx 78.69475\) when the simulation runs from 0 to 2000 s.

The comparison of various order moments along with the number density function obtained numerically is conducted with analytical results in Fig. 4. Analogous to the Section 4.1, the zeroth- and first-order moments predicted by the proposed scheme shows excellent agreement with the analytical results, that is, the number preservation and mass conservation properties hold for the proposed scheme as expected. The second-order moment calculated by the proposed scheme is also matching well with the analytical result (see Fig. 4b). Additionally, the number density function is plotted against the index of the cell reveals that the numerical results overlap with the analytical result (see Fig. 4c). To quantify the difference between the analytical and numerical values of the number density function, the weighted sectional errors are calculated and listed in Table 3. It shows similar results to the previous case as the errors decrease to approximately 50% when the proposed scheme is implemented on a refined grid of 60 nonuniform cells.

Fig. 4
figure 4

Comparison of various order moments and number density for binary breakage kernel, i.e., b(m,m) = 2/m with quadratic selection function S(m) = m2 for monodisperse initial condition

Table 3 Weighted error of number density for binary breakage kernel with quadratic selection function

In addition, the experimental order of convergence calculated for this case using uniform and nonuniform grids along with the L1 norm errors are listed in Table 4. It can be observed the prediction are very similar to Section 4.1 as the proposed scheme is exhibiting second-order convergence on both uniform and nonuniform grids. Moreover, it can also be observed that the L1 norm errors computed using uniform cells are very large compared to the nonuniform grids. However, these errors are reducing very fast when a refined grid is used. But, still, the proposed scheme with nonuniform grids perform better than the uniform grids (confirming the results provided by [9]).

Table 4 Experimental order of convergence for case II corresponding to monodisperse initial condition

4.3 Test case III

This section is devoted to validate the accuracy of the proposed scheme by taking into consideration a complex physically relevant breakage kernel of the following form:

$$ b(y,\xi) = \frac{12y}{\xi^{2}}\left( 1-\frac{y}{\xi}\right). $$

The quadratic selection function S(y) = y2 is considered for this case. The analytical solutions for moments and number density function for this kernel are not available in the literature. The accuracy for this kernel is measured by examining the experimental order of convergence computed using the expression given in [19]. The computation are run till time t = 200 s and attained the degree of breakage \(\displaystyle \frac {\mu (t)}{\mu (0)}\approx 12\). The computational size is considered as same as Section 4.1. It can be seen from Table 5 that the experimental order of convergence for the proposed scheme is similar to the previous two cases as it shows second-order convergence irrespective of the type of grid used to approximate these results.

Table 5 Experimental order of convergence for case III corresponding to monodisperse initial condition

5 Conclusions

In this work, a finite volume scheme has been proposed for approximating a fragmentation equation based on the idea of overlapping of the cells. The mathematical formulation of the proposed scheme is very simple and robust to implement on both uniform and nonuniform grids. The accuracy of the scheme has been verified for analytically tractable and physically tractable kernels. The qualitative and quantitative comparison shows that the scheme not only computed the zeroth- and first-order moments with high precision but also computed the second-order moment very accurately along with the number density function. The mathematical analysis of the proposed scheme has been also discussed in detail. The experimental order of convergence is compared with the theoretical observations in order to confirm that the proposed scheme is second-order convergent irrespective of the grid chosen for discretizing the given domain.