Convergence analysis of a finite volume scheme for solving non-linear aggregation-breakage population balance equations

This paper presents stability and convergence analysis of a finite volume scheme (FVS) for solving aggregation, breakage and the combined processes by showing Lipschitz continuity of the numerical fluxes. It is shown that the FVS is second order convergent independently of the meshes for pure breakage problem while for pure aggregation and coupled equations, it shows second order convergent on uniform and non-uniform smooth meshes. Furthermore, it gives only first order convergence on non-uniform grids. The mathematical results of convergence analysis are also demonstrated numerically for several test problems.


Introduction
The aggregation-breakage population balance equations (PBEs) are the models for the growth of particles by combined effect of aggregation and breakage. Each particle is identified here by its size, i.e. volume or mass. The equations we consider in this paper describe the time evolution of the particle size distribution (PSD) under the simultaneous effect of binary aggregation and multiple breakage. In binary aggregation, two particles combine together to form a bigger one whereas in breakage process, a big particle breaks into two or many fragments. There are many engineering applications, including aerosol physics, high shear granulation, highly demanding nano-particles and pharmaceutical industries etc., see Sommer et al. [28], Gokhale et al. [6] and references therein. Binary breakage is not sufficient for some of these applications, therefore, multiple fragmentation is preferred. The temporal change of the particle number density, f (t, x) ≥ 0, of particles of volume x ∈ R >0 at time t ∈ R >0 in a spatially homogeneous physical system undergoing an aggregation-breakage process is described by the following well known PBEs, see [23,32] ∂f (t, x) ∂t = 1 2 with initial data The first two terms on the right-hand side (rhs) are due to aggregation while the third and fourth terms model the breakage process. The two positive terms describe the creation of particles of size x and are called the birth terms for aggregation respectively breakage. The two negative terms describe the disappearance of particles of size x and are commonly called the death terms. The aggregation kernel β(x, y) ≥ 0 characterizes the rate at which two particles of volumes x and y combine together. It also satisfies the symmetry condition β(x, y) = β(y, x). The selection function S(ǫ) describes the rate at which particles of size ǫ are selected to break. The breakage function b(x, ǫ) for a given ǫ > 0 gives the size distribution of particle sizes x ∈ [0, ǫ[ resulting from the breakage of a particle of size ǫ. For the particular case of b(x, ǫ) = 2/ǫ, the multiple breakage PBE turns into the binary breakage PBE. The breakage function has the following important properties The functionN (x), which may be infinite, denotes the number of fragments obtained from the breakage of particle of size x. The second integral shows that the total mass created from the breakage of a particle of size x is again x. In aggregation-breakage processes the total number of particles varies in time while the total mass of particles remains conserved. In terms of f , the total number of particles and the total mass of particles at time t ≥ 0 are respectively given by It is easy to show that the total number of particles M 0 (t) decreases by aggregation and increases by breakage processes while the total mass M 1 (t) does not vary during these events. For the total mass conservation holds. However, for some special cases of β when it is sufficiently large compared to the selection function S, a phenomenon called gelation occurs. In this case the total mass of particles is not conserved, see Escobedo et al. [4] and further citations for details.
Mathematical results on existence and uniqueness of solutions of the equation (1) and further citations can be found in McLaughlin et al. [22] and W. Lamb [14] for rather general aggregation kernels, breakage and selection functions. In our analysis we consider them to be twice continuously differentiable functions. The PBEs (1) can only be solved analytically for a limited number of simplified problems, see Ziff [32], Dubovskii et al. [3] and the references therein. This certainly leads to the necessity of using numerical methods for solving general PBEs. Several numerical methods have been introduced to solve the PBEs. Stochastic methods (Monte-Carlo) have been developed, see Lee and Matsoukas [15] for solving equations of aggregation with binary breakage. Finite element techniques can be found in Mahoney and Ramkrishna [19] and the references therein for the equations of simultaneous aggregation, growth and nucleation. Some other numerical techniques are available in the literature such as the method of successive approximations by D. Ramkrishna [26], method of moments [18,21], finite volume methods [24,10] and sectional methods [8,12,30] to solve such PBEs.
A completely different numerical approach was proposed by Filbet and Laurençot [5] for solving aggregation PBEs by discretizing a well known mass balance formulation. They thereby introduced an application of the FVS to solve the aggregation problem. Further, Bourgade and Filbet [1] have extended their scheme to solve the case of binary aggregation and binary breakage PBEs and gave a convergence proof of approximate solutions in the space L ∞ (0, T ; L 1 (0, R)). For a special case of a uniform mesh they have shown error estimates of first order. The scheme has also been extended to two-dimensional aggregation problems by Qamar and Warnecke [25]. Finally it has been observed that the FVS is a good alternative to the methods mentioned above for solving the PBEs due to its automatic mass conservation property.
Since Bourgade and Filbet have considered aggregation with binary breakage problems on uniform meshes only. The objective here is to analyze such a FVS to solve the aggregation with multiple breakage PBEs on general meshes. We also demonstrate mathematically the missing stability and the convergence analysis of the FVS for simultaneous aggregation-breakage PBEs by following Hundsdorfer and Verwer [7] and Linz [17]. The mathematical results are verified numerically for several test problems on four different types of uniform and non-uniform grids.
This paper is organized as follows. First, we derive the FVS to solve aggregation-breakage PBEs. Then in Section 3 some useful definitions and theorems are reviewed from [7,17] which are used in further analysis of the method. Here we also discuss the consistency and prove the Lipschitz continuity of the numerical fluxes to get the convergence results. Later on the convergence analysis is numerically tested for several problems in Section 4. Further, Section 5 summarizes some conclusions. At the end of the paper one Appendix is provided which gives a bound on total number of particles for the aggregation-breakage terms.

Finite volume scheme
In this section a FVS for solving aggregation-breakage PBEs is discussed. Following Filbet and Laurençot [5] for aggregation, a new form of the breakage PBE is presented in order to apply the FVS efficiently. Then stability and convergence analysis will be discussed for the method.

Aggregation-breakage PBE in a conservative form
Writing the aggregation and breakage terms in divergence form enable us to get a precise amount of mass dissipation or conservation. It can be written in a conservative form of mass density The abbreviations agg and brk are used for aggregation and breakage terms respectively. The flux functions F agg and F brk are given by It should be noted that both forms of aggregation-breakage PBEs (1) and (4) are interchangeable by using the Leibniz integration rule. The concept of this conservative formulation of the PBE has been used in Tanaka et al. [29] and Makino et al. [20]. It should also be mentioned that the equation (4) reduces into the case of pure aggregation or pure breakage process when F brk (t, x) or F agg (t, x) is zero, respectively.
In the PBE (4) the volume variable x ranges from 0 to ∞. In order to apply a numerical scheme for the solution of the equation a first step is to fix a finite computational domain Ω :=]0, x max ] for an 0 < x max < ∞. Hence, for x ∈ Ω and time t ∈ (0, T ] where T < ∞, the aggregation and the breakage fluxes for the truncated conservation law for n, i.e. for are given as Here the variable n(t, x) denotes the solution to the truncated equation. We are given with initial data For further analysis, all the kinetic parameters β, S and b are considered to be two times continuously differentiable function, i.e.
From (11), there exists some non-negative constants Q and Q 1 depending on x max such that β(x, y) ≤ Q and b(x, y)S(y) ≤ Q 1 for x, y ∈]0, x max ].
Remark 2.1. The formulation we use here is a non-conservative truncation for the pure aggregation operator as F agg (t, x max ) ≥ 0 while it is mass conserving for the pure breakage equation, i.e. F brk (t, x max ) = 0. Hence, the combined formulation (7) is a non-conservative truncation as used by Bourgade and Filbet [1]. One could make a conservative truncation by replacing x max by x max − u in (8). This would give F agg (t, x max ) = 0. But it describes an artificial interruption of the aggregation process without a real physical justification. With our truncation particles that are too large leave the system.

Numerical discretization
Finite volume methods are a class of discretization schemes used to solve mainly conservation laws, see LeVeque [16]. For a semi-discrete scheme, the interval ]0, x max ] is discretized into small cells .., I, with where ∆x is the maximum mesh size. The representative of each size, usually the center of each cell x i = (x i−1/2 + x i+1/2 )/2, is called pivot or grid point. The FVS has been carried over to the discretization of such equations by instead of interpretingn i (t) as an approximation to a point value at a grid point, i.e. n(t, x i ), rather taking an approximation of the cell average of the solution on cell i at time tn Integrating the conservation law on a cell in space Λ i , the FVS is given as [16] x The term J − i+1/2 is called the numerical flux which is an appropriate approximation of the truncated continuous flux function F agg and/or F brk depending upon the processes under consideration. In case of a breakage process, the numerical flux may be approximated from the mass flux F brk as follows Using our assumptions that S ∈ C 2 (]0, x max ]), b ∈ C 2 (]0, x max ]×]0, x max ]) and applying the mid point rule we can rewrite (15) as Similarly for the aggregation problem, From Filbet and Laurençot [5], the above equation can be written as Here, the parameter I denotes the number of cells. The integer α i,k corresponds to the index of each cell such that Applying mid point approximation for the first term and Taylor series expansion of the second term about the point Let us denote the vector n := [n 1 , . . . , n I ] obtained by L 2 projection of the exact solution n into the space of step functions constant on each cell. It is worth to mention that this projection error can easily be shown of second order, see remark 3.3.3 in [11]. We also define the vectors ∆J agg (n) := [∆J agg 1 (n), . . . , ∆J agg I (n)] and ∆J brk (n) := [∆J brk 1 (n), . . . , ∆J brk Substituting the values of J agg i+1/2 and J brk i+1/2 from equations (19) and (16), respectively to get and By denoting the vectorn := [n 1 , . . . ,n I ] for the numerical approximations of the average values of n(t, x), the equation (14) can be rewritten as In order to retain the overall high accuracy, the semi-discrete scheme (23) can be combined with any higher order time integration method. It is worth to mention here that dealing with the pure cases of aggregation or breakage is easy by setting one of the two numerical fluxes is zero.

Convergence analysis
Before discussing the convergence of the semi-discrete scheme, let us review some useful definitions and theorems from [7,17] that will be used in the subsequent analysis. Let · denote the discrete L 1 norm on R I that is defined as In this work, we deal with this norm by interpreting the discrete data as step functions.
Definition 3.1. The spatial truncation error is defined by the residual left by substituting the exact solution n(t) = [n 1 (t), . . . , n I (t)] into equation (23) as The scheme (23) is called consistent of order p if, for ∆x → 0, It is important that our numerical solution remains non-negative for all times. This is guaranteed by the next well known theorem where we haveM ≥ 0 for a vectorM ∈ R I iff all its components are non-negative. Then the solution of the semi-discrete system (14) is non-negative if and only if for any vector n ∈ R I and all i = 1, . . . , I and t ≥ 0, Now we state a useful theorem from Linz [17] which we use to show that the FVS is convergent.
Theorem 3.4. Let us assume that a Lipschitz condition on J(n) is satisfied for 0 ≤ t ≤ T and for all n,n ∈ R I where n andn are the projected exact and numerical solutions defined in (7) and (23), respectively. More precisely there exists a Lipschitz constant L < ∞ such that holds. Then a consistent discretization method is also convergent and the convergence is of the same order as the consistency.
Proof. A more general result is proven in Linz [17].
Due to Theorem 3.4, for the convergence of our scheme it remains to show that the method is consistent and the Lipschitz condition (26) is satisfied by the fluxes.

Consistency
The following lemma gives the consistency order of the FVS for aggregation-breakage PBEs.
Lemma 3.5. Consider the function S ∈ C 2 (]0, x max ]) and b, β ∈ C 2 (]0, x max ]×]0, x max ]). Then, for any family of meshes, the consistency of the semi-discrete scheme (23) is of second order for the pure breakage process, i.e. with ∆J agg (n) = 0. For the aggregation and coupled processes, the scheme is second order consistent on uniform and non-uniform smooth meshes while on oscillatory and random meshes it is first order consistent.
Proof. The spatial truncation error (25) is given by Integrating (7) over Λ i and applying the mid-point rule in the time derivative term, we interpret Substituting this into the equation (27) and using (20) give the following form Let us now begin with We now use Taylor series expansion of the functions K x i±1/2 (ǫ) := n(t, ǫ) x i±1/2 0 ub(u, ǫ) du about x k and further rearrangement of terms yield σ brk i (t) as Applying the mid-point rule, it should be noted that Thus we obtain σ brk i (t) = O(∆x 2 ). Hence, for the pure breakage process, the consistency of the semi-discrete scheme (23) is two which is determined by using (24) as independently of the type of meshes.
Due to the non-linearity of the aggregation problem, it is not easy to determine the consistency order on general meshes and therefore, we evaluate it on various meshes separately. The results can be combined to the results of breakage process to give the consistency of the coupled processes. We know from (17) Applying the mid-point rule, it should again be noted that Therefore, by defining LHS : Substituting the values of L x i±1/2 (x j ) yield (leaving the third order terms) Now, I 1 is equivalent to Applying the mid-point approximation for the second term, we figure out Similarly, we estimate Subtracting the third term from I 2 to I 1 gives By using Lemma 3.6 which is stated in the next section, the summation over k is finite in this term. Hence, the rhs of this equation becomes of order O(∆x 3 ) and can be omitted. Therefore, Open the Taylor series about the points x α i,j −1 in I 3 and x α i−1,j −1 in I 4 as well as by using the relation (19), we finally obtain Now the consistency order on four different types of meshes are evaluated:

Uniform mesh
Let us assume that the first mesh is uniform, i.e. ∆x i = ∆x for all i. In this case x i+1/2 − x j and x α i,j −1 become the same and are equal to the pivot point x i−j+1 . Similarly, Applying the Taylor series expansion of the function ) about the point x α i−1,j −1 in the first term on the rhs of the equation (30) to get uniform mesh non-uniform mesh for example x = exp(ξ) Figure 1: Non-uniform smooth mesh.
Further by facilitating the integrals and using the relation (31), we have Hence, σ agg i (t) = O(∆x 2 ) and so the order of consistency is given by using (24) as Therefore, the scheme is second order consistent on uniform grids.

Non-uniform smooth mesh
A smooth transformation from uniform grids leads to such meshes. In this case grids are assumed to be smooth in the sense that ∆x where ∆x is the maximum mesh width. For example, let us consider a variable ξ with uniform mesh and a smooth transformation x = g(ξ) to get non-uniform smooth mesh, see Figure 1.
For the analysis here, we have considered the exponential transformation as x = exp(ξ). Such a mesh is also known as a geometric mesh, i.e. x i+1/2 = rx i−1/2 with r = exp(h). The termh is the width of the uniform grid. Here again we achieve second order consistency.
Equation (30) can be rewritten by setting j = j − 1 in second term as Further it can be rewritten as It can further be simplified as Further notice that x i+1/2 − x j = r(x i−1/2 −x j−1 ) and therefore α 1 = α 2 +1. Again by using the condition ∆x j −∆x j−1 = O(∆x 2 ), we determine ∆x 2 α 1 − ∆x 2 α 2 = O(∆x 3 ). Now, to get a second order consistency of the scheme, it is remained to show that or equivalently, Let us consider ξ 1 , ξ 2 are corresponding points in the uniform mesh for x α 2 and x i−1/2 − x j−1 , respectively. Consider h 1 = ξ 2 − ξ 1 which is given as Similarly, taking h 2 = ξ 4 − ξ 3 where ξ 3 and ξ 4 are the points in the uniform mesh corresponding to the points x α 1 and x i+1/2 − x j , respectively, we evaluate Finally, the equation (32) can be estimated by using Taylor series expansion as Hence, by using (28) and (24) the order of consistency for the pure aggregation process is two for the smooth meshes x i+1/2 = rx i−1/2 .

Oscillatory and random meshes
A mesh is known to be an oscillatory mesh, if for r > 0(r = 1) it is given as From the equation (30), it is clear that the first two terms on the rhs can not be cancel out for an oscillatory or a random mesh. Therefore, σ agg i (t) = O(∆x) and so the accuracy of the semi discrete scheme (23) is one by using the relation (24) on such meshes. Now for the coupled aggregation and breakage problems, the local truncation error of each process can be combined and give second order consistency on uniform and non-uniform smooth meshes whereas it is of first order on the other two types of grids.

Lipschitz continuity of the fluxes
To prove the Lipschitz continuity of the numerical flux J(n) in (23), the following three lemmas are used. Lemma 3.6. Let us assume that the points x i+j− 1 2 − x k for given i, k and j = 1, 2, . . . , p where p ≥ 2 lie in the same cell Λ α for some index α. We also assume that our grid satisfies the quasi-uniformity condition for some constant C(independent of the mesh size). Then p is bounded by C + 1.
In the next two lemmas the boundedness of the total number of particles for the aggregation and multiple breakage equations are discussed.
Lemma 3.7. Let us assume that the kernels β, S and b satisfy the boundedness condition (12). Then the total number of particles for the continuous aggregation-breakage equation (7) is bounded by a constant C T,xmax > 0 depending on T and x max , namely Proof. The proof can be found in Appendix A.
Lemma 3.8. Under the same assumptions on β, S and b considered in the previous lemma, we have boundedness of the total number of particles for the discrete aggregation-breakage equation (14) by using the finite volume scheme. The bound in this case is again C T,xmax as before, i.e.
provided that the initial dataN (0) and N (0) are the same.
Proof. The proof has been given in Appendix A. Now, the Lipschitz continuity of the numerical flux J(n) defined as in (23) is shown.
Lemma 3.9. Let us assume that our grid satisfies the quasi-uniformity condition (34). We also assume that the kernels β, S and b satisfy the bounds (12) which are β ≤ Q and bS ≤ Q 1 . Then there exists a Lipschitz constant L := (4C + 6)QC T,xmax + 2Q 1 x max < ∞ for some constants C, C T,xmax > 0 such that holds.
Proof. From (23), we have the following discretized form of the equation To prove the Lipschitz conditions on J(n), it is sufficient to find the Lipschitz conditions on ∆J agg (n) and ∆J brk (n) separately. For the aggregation, Substituting the value of ∆J agg i (n) from the equation (21) yields Now the terms S i , i = 1, · · · , 4 in (38) are evaluated one by one. First the term S 1 is simplified which may be estimated Since k < i implies that x k < x i . Using the relation xy −xŷ = 1/2[(x−x)(y +ŷ)+(x+x)(y −ŷ)], bound β(x, y) ≤ Q and setting N i = n i ∆x i give Open the summation for each i, we obtain Having Lemmas 3.7 and 3.8, which say that the total number of particles is bounded by a constant C T,xmax , S 1 is further simplified as S 1 ≤ 2QC T,xmax n −n .
Now the term S 2 is calculated from (38) which is taken as Further simplifications as in the previous case yield Changing the order of summation gives By using the Lemma 3.6 which shows that the number of repetition of index in a cell is finite and bounded by some constant C, we obtain S 2 ≤ 2CQC T,xmax n −n . The same bound on S 3 is achieved because the only difference is that the index i − 1 is used instead of i.
Finally the expression S 4 from (38) can be written as Further simplification gives S 4 ≤ 4QC T,xmax n −n . Adding all the results from S 1 , S 2 , S 3 and S 4 yields with a Lipschitz constant L 1 = (4C + 6)QC T,xmax .
Similarly, for the breakage problem, we have By using the equation (22), the above equation reduces to Since x j < x i for j < i and having bS ≤ Q 1 from (12), the above can be simplified as Therefore, the following is obtained with a Lipschitz constant L 2 = 2Q 1 x max . Hence, the Lipschitz conditions for J(n) with a Lipschitz constant L = (4C + 6)QC T,xmax + 2Q 1 x max is shown.
Hence, by Theorem 3.4 the order of convergence of the FVS for the aggregation or breakage or coupled processes is same as the order of consistency which we have seen before in Lemma 3.5.

Numerical Results
The mathematical results on convergence analysis are verified numerically for pure aggregation, breakage and also for the combined processes considering several test problems. All numerical simulations below were carried out to investigate the experimental order of convergence (EOC) on four different types of meshes discussed in the next subsection.
If the problem has analytical solutions, the following formula is used to calculate the EOC EOC = ln(E I /E 2I )/ ln (2).
Here E I and E 2I are the discrete relative error norms calculated by dividing the error N −N by N where N,N are the number of particles obtained mathematically and numerically, respectively. The symbols I and 2I correspond to the number of degrees of freedom. Now, in case of unavailability of the analytical solutions, the EOC can be computed as whereN I is obtained by the numerical scheme using a mesh with I degrees of freedom.
Before going into the details of the test cases, in the following subsection we discuss briefly four different types of uniform and non-uniform meshes where global truncation errors are obtained numerically. These meshes have also been used in J. Kumar and Warnecke [9].

Meshes
Uniform mesh: A uniform mesh is obtained when ∆x i = ∆x for all i.
Non-uniform smooth mesh: We are familiar with such a mesh from the previous section and Figure 1. For the numerical computations, a geometric mesh is considered.
Oscillatory mesh: The numerical verification has been done on an oscillatory mesh by taking r = 2 in the equation (33). In this case, the EOC is evaluated numerically by dividing the computation domain into 30 uniform mesh points initially. Then each cell is divided by a 1:2 ratio on further levels of computation.
Random mesh: Similar to the previous case, we started again with a geometric mesh with 30 grid points but then each cell is divided into two parts of random width in the further refined levels of computation. Here, we performed ten runs on different random grids and the relative errors are measured. The average of these errors over ten runs is used to calculate the EOC.

Pure aggregation
Test case 1: The numerical verification of the EOC of the FVS for aggregation is discussed by taking two problems, namely the case of sum and product aggregation kernels. The analytical solutions for both problems taking the negative exponential n(0, x) = exp(−αx) as initial condition has been given in Scott [27]. Hence, the EOC is computed by using the relation (41). Table 1 shows that the EOC is 2 on uniform and non-uniform smooth meshes and is 1 on oscillatory and random grids in both cases. The computational domain in this case is taken as [1E − 6, 1000] which corresponds to the ξ domain [ln(1E − 6), ln(1000)] for the exponential transformation x = exp(ξ) for the geometric mesh. The parameter α = 10 was taken in the initial condition. The simulation result is presented at time t = 0.5 and t = 0.3 respectively for the sum and the product aggregation kernels corresponding to the aggregation extentN (t)/N (0) ≈ 0.80.

Pure breakage
Test case 2: Here, the EOC is calculated for the binary breakage b(x, y) = 2/y together with the linear and quadratic selection functions, i.e. S(x) = x and S(x) = x 2 . The analytical solutions for such problems have been given in Ziff and McGrady [31] for a mono-disperse initial condition of size unity, i.e. n(0, x) = δ(x − 1). Hence, by using the formula (41), we observe from the Table 2 that the FVS is second order convergent on all the grids. The computational domain in this case is taken as [1E − 3, 1]. Since the rate of breaking particles taking quadratic selection function is less than that of linear selection function, we take t = 100, 200 for linear and quadratic selection functions, respectively. The time has been chosen differently for both the selection functions to have the same extent of breakageN (t)/N (0) ≈ 22.

Test case 3:
Now the case of multiple breakage with the quadratic selection function S(x) = x 2 is considered where an analytical solution is not known. Therefore, the EOC is calculated using (42). For the numerical simulations, the following normal distribution as an initial condition is taken The computations are made for two breakage functions considered by Diemer and Olson [2] and Ziff [32], respectively In case(i) the relation y 0 b(x, y)dx = p holds where p gives the total number of fragments per breakage event. The parameter c ≥ 0 is responsible for the shape of the daughter particle distribution, see also [28]. The numerical solutions are obtained using p = 4, c = 2. The second breakage function gives ternary breakage. For the numerical simulation the minimum and maximum values of x are taken as 1E − 3 and 1 respectively. The time t = 100 is set to get the breakage extentN (t)/N (0) ≈ 22 in case(i) while t = 150 is used for case(ii). As expected from the mathematical analysis, we again observe from the Table 3 that the FVS shows convergence of second order on all the meshes. The computations for higher values of p up to 19 are also tested and observed that there is no marked difference in the EOC.

Coupled aggregation-breakage
Test case 4:   Finally, the EOC is evaluated for the simultaneous aggregation-breakage problem considering a constant aggregation kernel β(x, y) = β 0 and breakage kinetics b(x, y) = 2/y, S(x) = x. The analytical solutions for this problem are given by Lage [13] for the following two different initial conditions This is a special case where the number of particles stays constant. The later initial condition is a steady state solution. For the simulation the computational domain [1E − 2, 10] with N 0 = x 0 = 1 and time t = 0.3 is taken. From Table 4, we find that the FVS is second order convergent on uniform and non-uniform smooth meshes and it gives first order on oscillatory and random meshes using (41). It should be mentioned that the computation has also been done for the product aggregation kernel β(x, y) = xy and the linear selection function S(x) = x taken together with two different general breakage functions as stated in the previous section. Analytical solutions are not available for such problems and so the EOC was calculated using (42). We observed again that the FVS shows similar results of convergence for these meshes.

Conclusions
In this article the convergence analysis of the finite volume techniques was studied for the non-linear aggregation and multiple breakage equations. We showed the consistency and then proved the Lipschitz continuity of the numerical fluxes to complete the convergence results. This investigation was based on the basic existing theorems and definitions from the book of Hundsdorfer and Verwer [7] and the paper of Linz [17]. It was noticed that the scheme was second order convergent for a family of meshes for the pure breakage problem. For the aggregation and combined processes, it was not straightforward to evaluate the consistency and the convergence error on general meshes. This depended upon the type of grids chosen for the computations. Moreover, in these cases the method gave second order convergence on uniform and non-uniform smooth meshes while on non-uniform grids it showed only first order. The mathematical results of convergence analysis were verified numerically on several meshes by taking various examples of pure aggregation, pure breakage and the combined problems.

A Bound on total number of particles
We give the proof of Lemmas 3.7 and 3.8 in Appendices A.1 and A.2, respectively.

A.1 Continuous aggregation and multiple breakage equation
Proof.
[ Lemma 3.7] Integrating the equation (7) with respect to x from 0 to x max gives From the equations (8) and (9), we know that ∂ ∂x (F agg (t, x)) = ∂ ∂x Since x ≥ u for the first integral, this implies that u/x ≤ 1. Substituting x = z + u such that dx = dz, the above can be rewritten as Notice that the first two integrals combined give a negative value. Using the relation (3) of the breakage function in the last integral and due to negativity From the bounds (12) we know that bS ≤ Q 1 . Estimating v ≤ x max leads to Therefore, the total number of particles is bounded and the bound is given as

A.2 Discrete aggregation and multiple breakage equation
Proof.
[ Lemma 3.8] Multiplying the equation (14) by ∆x i /x i and summing with respect to i gives We write out the summation over i of the aggregation fluxes J agg i±1/2 to get For the breakage fluxes J brk i±1/2 in (49) we substitute the definition (16). Introducing the notationsN i (t) = n i (t)∆x i andN (t) = I i=1N i (t) ensure Due to positivity of J agg i+1/2 for all i and J agg 1/2 = 0, we estimate Changing the order of summation for the first term and the summation indices in the second term yield Since i < k implies that 1 − x i /x k < 1. Having the bound bS ≤ Q 1 gives dN (t)/dt ≤ x max Q 1N (t). Therefore, the following bound is obtained on the total number of particles by using the FVS aŝ which is the same bound as explained in the previous lemma, providedN (0) = N (0).