Rate of convergence and stability analysis of a modified fixed pivot technique for a fragmentation equation

This study presents the convergence and stability analysis of a recently developed fixed pivot technique for fragmentation equations (Liao et al. in Int J Numer Methods Fluids 87(4):202–215, 2018). The approach is based on preserving two integral moments of the distribution, namely (a) the zeroth-order moment, which defines the number of particles, and (b) the first-order moment, which describes the total mass in the system. The present methodology differs mathematically in a way that it delivers the total breakage rate between a mother and a daughter particle immediately, whereas existing numerical techniques provide the partial breakup rate of a mother and daughter particle. This affects the computational efficiency and makes the current model reliable for CFD simulations. The consistency and unconditional second-order convergence of the method are proved. This demonstrates efficiency of the method over the fixed pivot technique (Kumar and Warnecke in Numer Math 110(4):539–559, 2008) and the cell average technique (Kumar and Warnecke in Numer Math 111(1):81–108, 2008). Numerical results are compared against the cell average technique and the experimental order of convergence is calculated to confirm the theoretical order of convergence.


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.

Fragmentation equation
The temporal change of particle number density, n(t, v) ≥ 0 of volume v ≥ 0 at time 0 ≤ t ≤ T < ∞ in a spatially homogeneous physical system undergoing fragmentation is written as (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 ζ(y) produced during the breakup of particles size y. The last expression, indicates the symmetric nature of the breakage function.
Additionally, p (v, y) is the partial breakage kernel with an incorporated binary size distribution, which is mirror symmetrical about v = y 2 and defined by the product of the total breakup rate and the daughter size distribution p (v, y) = t (y) β(v, y) and t (y) = 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 (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.

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 where θ is the step function θ(v) := 0, when v is false, 1, when v is true. Thus the first term of the main Eq. (1.1) is redefined as . (1.6) To summarize, the first component B 1 (v) represents the accumulation of particle sizes smaller than y 2 , and the second component B 2 (v) represents the accumulation of particle sizes greater than y 2 , through the complementary particles x.
Mass conservation: We first focus on the term B 1 (v). Changing the order of integration Now θ v ≤ y 2 represents the accumulation of particle sizes smaller than y 2 , so we can write We next consider B 2 (v), and note that it represents the accumulation of particles v which are larger than y 2 . Therefore, x = y − v ≤ y 2 , and thus θ x ≤ y 2 = 1. In the first step, we change the order of integration Therefore the mass conservation law follows. Similarly we can calculate the zeroth moment as

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.

Computational domain and numerical scheme
The midpoint of each cell i is denoted by x i , and the cell width is denoted by Accordingly, the discrete number of particles in each cell (N i ) calculated from the number density function reads as To derive the mathematical formulation of the scheme, the birth term is integrated over the size range Applying several approximations, the birth term of the modified fixed pivot technique is written aŝ Here,N i is the numerical approximation of the number of particles in each cell, ν i ( j), Y i jk are the birth-modification factors which assign the particles to their nearest representatives, and are defined as Similarly, integrating the death term over each subinterval giveŝ Note that Liao et al. [11] have written the discrete partial breakage kernel p as follows: Under the above consideration, the modified fixed pivot technique of [11] is written as

Zeroth order moment consistency
During fragmentation of x j , the distribution of daughter particles is based upon their sizes relative to 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 x j 2 and get The distribution of particles larger than x j 2 is achived through the complementary particles x ≤ x j 2 . Due to symmetry of the breakage function, x falls in the same interval 1 ≤ j ≤ j, and thus we get Combining the two cases, we get Hence, the time evolution of the discrete zeroth moment is calculated as

Mass conservation law
First consider the distribution of the particles of size smaller than x j 2 .
Now the argument of cell allocation for particles follows similarly to the previous description, that is, consider that x j 2 falls in some cell j such that 1 ≤ j ≤ j. Therefore, For the distribution of particles larger than x j 2 , we compute Combining the two cases, we get Hence, time evolution of the discrete first moment is calculated as
Then the semi-discrete system (2.3) is expressed as where A is the I × I upper triangular matrix Throughout the subsequent discussion, we work with the discrete L 1 −norm defined as follows 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 N = {N 1 , N 2 , . . . , N I } into the discrete system as The scheme is called consistent of order q if, for x → 0, Definition 3.2 Consider a matrix A := a i j ∈ R m×m . Then for the L p vector norm, the corresponding logarithmic norm of A is defined as Moreover, the logarithmic norm of a matrix A corresponding to the L 1 norm on R m is given byμ holds for all grids with constant K ≥ 1 and ω ∈ R, both independent of x. Proof The proof of the above theorem can be found in [4] (Theorem 7.2, p. 117).
Proof From the definitions of all the functions and the matrix A , it can be easily seen that α i j ≥ 0 for all i = j. Thus recalling Theorem 3.2, we deduce that the solution of the semi-discrete system (3.1) is nonnegative.

Consistency
Let C 2 [a, b] denote the space of twice continuously differentiable functions on the nonempty bounded closed interval [a, b].
. 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 Integrating the birth term B (v) over each cell i we get that Applying quadrature rules for the outer integrals on the right-hand side, we get Again we apply the midpoint rule to the outer integrals of the first two terms and get 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 Therefore, Here, By the application of the midpoint rule, the definition of ν i ( j), Y i jk and other simple computations, we get We now proceed to estimate the error terms E 1 , E 2 and E 3 . In general, the definition of β yields that β(x, y) = 0 whenever x ≥ y.
Therefore, the definition of E 1 gives

Now for estimating E 2 , limits of the integrals give
Hence, we get For the death term, direct application of midpoint quadrature formula leads to Hence, This proves that the scheme is second order consistent irrespective of the choice of the spatial mesh.

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 Since all non-diagonal entries of the matrix A are nonnegative, the above logarithmic norm takes the following form Suppose that x j 2 falls in some interval j satisfying 1 ≤ j ≤ j. Then the definition of ν i ( j) redefines the formulation as Again from the definition, the lower boundary of Y i jk is given by x i−1 ≤ x j − x k ≤ x i which implies that the indices i and k are together replaced by a new index i which ranges from j + 1 → j. Therefore, further replacing i by i, we get Therefore the logarithmic normμ 1 (A) ≤ 2ω, and consequently by using Theorem 3.1 we get exp(t A) ≤ exp(2ωt) which ensures the stability of the scheme with stability constant 2ω.
Lipschitz condition: Consider the problem written as where the numerical flux is defined by for all 0 ≤ t ≤ T and N,N ∈ 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 Supposing that x j 2 falls on some interval 1 ≤ j ≤ j, we get where ω := max t (x). Therefore the Lipschitz constant ξ := 3ω < ∞.

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 (x) can be calculated using the zeroth and first moments from the following relation:x (4.1) The absolute maximum relative error in the moments is calculated using the following expression: where M i andM 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: Moreover σ i (t) provides the relative weighted sectional error in the number density function over the whole volume domain for i = 0.

Binary breakage kernel with linear selection function corresponding to exponential initial condition
For comparing the results, the binary breakage kernel with linear selection function ( t (x) = 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. 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  Table 2. Once again, the MFPT provides the numerical results with higher accuracy than the cell average technique.

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 ( t (x) = x 2 ) is chosen corresponding to the monodisperse initial condition n(0, x) = δ(x − 1). The computational grid and time required for running the simulation are taken to be same as in the previous case. 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.

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:   The subscript I corresponds to the degrees of freedom, and the relative L 1 error is computed by N −N N . The EOC is calculated for the binary breakage kernel with linear selection function corresponding to a monodisperse initial condition n(0, x) = δ(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].

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.
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/.