Finite volume approach for fragmentation equation and its mathematical analysis

This work is focused on developing a finite volume scheme for approximating a fragmentation equation. The mathematical analysis is discussed in detail by examining thoroughly the consistency and convergence of the numerical scheme. The idea of the proposed scheme is based on conserving the total mass and preserving the total number of particles in the system. The proposed scheme is free from the trait that the particles are concentrated at the representative of the cells. The verification of the scheme is done against the analytical solutions for several combinations of standard fragmentation kernel and selection functions. The numerical testing shows that the proposed scheme is highly accurate in predicting the number distribution function and various moments. The scheme has the tendency to capture the higher order moments even though no measure has been taken for their accuracy. It is also shown that the scheme is second-order convergent on both uniform and nonuniform grids. Experimental order of convergence is used to validate the theoretical observations of convergence.


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: subject to a initial condition 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: and ξ 0 yb(y, ξ )dy = ξ, ∀ ξ > 0. (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: 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 firstorder 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.

Numerical scheme
For developing the numerical method, the upper limit present in the first integral (2) must be replaced by a positive number say, y max . Thus, the equation takes the following form: Further, it is necessary to discretize the continuous computational domain from 0 to y max into L number of cells as shown in Fig. 1 where y L+1/2 ≤ y max . The lower and upper boundaries of the j th cell is denoted by y j −1/2 and y j +1/2 , respectively, for j ∈ N whose mean is y j = y j −1/2 + y j +1/2 2 . The step size of j th cell is given by y j = y j +1/2 − y j −1/2 . Now let us discretize the time domain as t n+1 = t n + t n for all n ∈ N. Further assume that the average value of f at any time t n in j th cell is the approximation of the function f y j , t n . Further, the function f can be discretized as f j (t) = f y j , t + O y 2 j 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 df j dt = B j − D j , for j = 1, 2, · · · , L.
subject to a new initial condition Here the birth and death terms are written as Consider the birth term (8) and simplify the equation by changing the order of integration, as shown below: The midpoint quadrature approximation is then employed to both integrals of (10) with respect to y, yielding: In a similar manner, the discretize form of the death term (9) can be simplified to Assume that in a j th cell,f j is the numerical approximation of the number density function f j . As a result, the discrete equations for approximating the number density function are as follows: Here where p j k = y k , when k = j, y k+1/2 , otherwise. Now we will define a proportionality constant in detail to introduce the notion of overlapping: For a uniform grid, the value of k j 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, k j can be calculated based on the fact that when a particle having properties, say y k of a kth cell breaks, it form two particles having properties y k−j and y j for a nonuniform grid. Then the lower and upper bounds of the newly formed particle after fragmentation mechanism becomes y j −1/2 = y k−1/2 − y (k−j)−1/2 and y j +1/2 = y k+1/2 −y (k−j)+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, y j −1/2 > y k−1/2 − y (k−j)−1/2 and y k+1/2 − y (k−j)+1/2 < y j +1/2 . • When the domain of newly born particle of size y k − y k−j will totally fall inside the cell, that is, The schematic representation of all possible cases of overlapping of the cells is shown in Fig. 2.
Mathematically, the proportionality constant for overlapping of the cells is calculated using the relation given below: where ∧ k j = min(y j +1/2 , y k+1/2 − y (k−j)+1/2 ) and ∨ k j = max(y j −1/2 , y k−1/2 − y (k−j)−1/2 ). Here, ∧ k j and ∨ k j define the bounds of the intersection of the cells k and (k − j) with cell j . The proportionality constant k j describes the extent of overlapping of the newly formed particle with cell j . The equality holds when the intersection is empty k j = 0 and when the newly born particles domain falls completely inside the cell j k j = 1 . Using the aforementioned relations and expressions, the formulation (12) becomes 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 Here the weights responsible for the preservation of total number of particles and conservation of the total mass in the system are given by and The values of weights ω b 1 and ω d 1 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 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: However, the numerical scheme holds the mass conservation property when it satisfies the following condition: Equations (20) and (21) can be easily obtained from (7) by multiplying with y j 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 Now changing the order of the summation, we have Substituting the values of weight from (18), the above equation gives Further replace the values of the weight from (17) and after simplification, it gives 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 y j on both sides and take summation over j will give The above can also be rewritten as Substituting the values of ω d k from (18) and after simplification, we get This further gives 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.

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 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 Here, the birth term (B), death term (D) and J ∈ R I are the functions off with the componentsB andD Therefore, the final form of equation is To conduct the convergence of the discrete system, first it is necessary to define the norm L 1 which is given by Now, let us give some important definitions which will be required in conducting the convergence analysis of the proposed finite volume scheme. f 2 (t), . . . , f L (t)] in the discrete system of equations, that is,

Definition 1 Spatial Truncation Error σ (t) is obtained by substituting the exact solution
The numerical scheme is said to be consistent of order p, if y → 0 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 (t) = f (t) −f (t). The numerical scheme is said to be convergent of order p if, for y → 0,

Proposition 3 Let us suppose that a Lipschitz condition on J (f ) is satisfied for
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.
such that the Lipschitz condition on J is satisfied for all f andf ∈ R L , that is, Proof Further consider the definition of the norm defined in (33) which further simplifies to Now simplify the first term by changing the order of the summation, we get After substituting the value of ω b k and using the fact that ω b 1 = 0 in the above equation, we have Rearranging the above equation, we get 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 Hence, Additionally, for k = 2, 3, · · · , L, we have (39) where M 2 is a constant, satisfying 0 < M 2 = min k∈{2,3,··· ,L} Using the above relations in (37), we get Here δy = min i y i . The above relation can further be simplified to where M 3 = 2 W max y∈{2,3,··· ,y max } [v(y)S(y)] ≤ ∞ and W = 2 1 + M 1 M 2 Ky max . Further, consider the second term T 2 from the (35), Using the fact that y j ≤ y k ∀ j = 1, 2, 3, · · · , k, Now changing the order of the integration leads to This gives Using the (43) and (44), the (35) takes the following form: where M = 2M 3 < ∞ 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. 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.
Consistency Using the Definition 1, the j th component of the spatial truncation error can be written as Further using (7) and (32), the above relation becomes Now consider the first term I 1 of the (46), Combining the terms, we get Further estimate the order of the term 1 − k j ω b k as follows: Using the mass conservation property of the breakage function defined in (3) gives 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 Using the above relation, the (48) implies Similar to the term I 1 , now let us discuss the order of consistency for the term I 2 of (46) Substituting the value of ω d j from (18) in the above equation, we get It can be noted that Using y j from the above equation in (51) will give As proved earlier that 1 − ω b j = O y 2 , we have Hence, from (46), (50) and (55), we get 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.

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: The comparison is enhanced by quantifying the weighted sectional errors exist in the number density function which can be estimated using the following relation: 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]: where E I and E 2I describe the L 1 error norm calculated by Here, f exc j and f num j describe the number of particles obtained analytically and numerically, respectively, and the symbols L and 2L correspond to the number of degrees of freedom.

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 y min = 10 −9 and ends with y max = 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, μ(t) μ(0) ≈ 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 zerothand 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.
In order to quantify the errors in the number density function, the weighted sectional errors calculated using the expression (57) 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 L 1 norm (59) decreases enormously as the number of cells increased in the domain.

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  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. In addition, the experimental order of convergence calculated for this case using uniform and nonuniform grids along with the L 1 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 L 1 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]).

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: The quadratic selection function S(y) = y 2 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 μ(t) μ(0) ≈ 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.

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.
Funding Open Access funding provided by the IReL Consortium.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.