Numerical solutions for multidimensional fragmentation problems using finite volume methods

We introduce a finite volume scheme for approximating a general multidimensional fragmentation problem. The scheme estimates several physically significant moment functions with good accuracy, and is robust with respect to use of different nonuniform daughter distribution functions. Moreover, it possess simple mathematical formulation for defining in higher dimensions. The efficiency of the scheme is validated over several test problems.


1.
Introduction. Fragmentation or breakage, comprising all processes in which two or more objects are created from one initial object, occurs in many natural and industrial processes, for example: cell division, bubble break-up, milling and grinding of particulate materials [9,19,22,33].
The fragmentation behaviour, i.e. when does fragmentation occur and to what extent, and the properties of the newly created fragments, is determined by a number of internal properties and external conditions. External conditions comprise for instance the medium in which the fragmentation process occurs, e.g. a dry or wet milieu, or external conditions that create forces on the initial objects, force or flow fields. Internal properties of an object that determine the fragmentation behaviour are for example: Object size, object shape and fractal dimension, internal porosity (void space), and, in case of heterogeneous materials, the composition. If the initial object is an agglomerate, i.e. consisting of a number of primary particles connected via attractive forces or solid or liquid bridges, then the fragmentation behaviour depends also on the number and individual strength of the contacts between the primary particles.
Fragmentation processes are therefore multidimensional, determined by geometric as well as material properties. Modeling of these processes, for instance for process design and optimisation, is inherently difficult, as they are also of a multiscale nature: fragmentation at single-object level and fragmentation at the level of the whole population of objects, for example all bubbles in a bubble column or particles in a grinding process. Due to its industrial importance, several studies have been performed to describe fragmentation at single-object level, for instance [7,8,12,30], using Discrete Element methods (DEM) or from experimental data [22,35,37,38]. Additionally, fragmentation has been studied at the level of populations, utilising the population balance approach (PBM) [2,14,25,40], allowing the estimation of product properties, e.g. size distribution of fragments, which can then be used in reverse to design process equipment and set-up operating conditions. However, due to the multidimensional character, obtaining numerical results is still a computational issue (as will be outlined in the following), complicating the design and operation of fragmentation processes.
In order to approach a solution by incorporating different particle properties, we formalize the problem under consideration: for any integer d ≥ 1, let R d + be the d−dimensional real space with elements x = (x 1 , x 2 , . . . , x d ) such that x r ≥ 0 for all r = 1, 2, . . . , d. In a general multidimensional system, the components of x represent particle properties like size, energy, moisture content, shape factors etc. In this context, there occur certain fragmentation events where one can describe the fragmentation model by a unique variable, say size. Under such scenario all the other variables are eliminated by some appropriate mathematical treatment. Thus reducing the system into its one-dimensional counterpart [3]. However, this elimination of variables is not possible for any multidimensional fragmentation. Therefore a density function, f ( x, t) is needed to describe the particle property distribution x at time t(≥ 0). Thus the multidimensional breakage population balance equation is written as which is supplemented by the initial data The functions S( y) and b( x| y) in the above equation, respectively denote the selection rate and distribution of fragmented daughter particles. In general, each of f , S and b are non-negative functions. Here, a function is considered to be non-negative whenever each component of that function is non-negative. The integral in (1) basically represents the following product of integrals over the internal coordinates, For the PBEs, different moments of the particle density function play significant roles as some of them correspond to important physical properties like total number, mass, energy, shape factor of the particles etc., [3,10,16,32]. Let us now formally introduce the moments of the number density function f ( x, t). Considering p := d r=1 p r , where each p r is a non-negative integer, the p−th order moment of f ( x, t)

81
is given by Similar to the one-dimensional problems, the zeroth moment M 0,...,0 (t) corresponds to the total number of particle present in the system and the temporal change of zeroth moment can be written as where ν( y) := y 0 b( x| y) d x, denotes the number of fragments produced during a breakage event.
Since in a one-dimensional model, the particle property is defined by a single variable namely, size. Therefore, the first moment corresponds to the total mass of the particulate system. Similarly, in a multidimensional system when each component of x represents size, say length, the d−th order moment M 1,...,1 (t) represents the hypervolume of the particles. This hypervolume is basically the geometry of the d−dimensional particles. For example, in two-dimensions if the two internal coordinates of x represents length and width of the particle, then hypervolume corresponds to the total area of the particle. A detailed description on the concept of hypervolume can be found in the articles of [1,26,32]. However the above description cannot be qualitatively inherited for a general multidimensional system. Therefore, when each component of x represents different particle properties, the first moment corresponding to a certain property defines the total content of that particular property in the system. Mathematically, the first moment M 0,...,1,...,0 (t), 1 in the r−th position, corresponds to the total content of the x r −property. From relation (3), we observe that there are d first-order moments. Therefore conservation of the total content of the particle properties requires the conservation of the first-order moments taken together [4].
In this regard, it is important to mention that the number of conserved moments depends upon the choice of the selection function S( y) and the breakage function b( x| y). In various industrial sectors an accurate estimation of the zeroth and the first moments is very important. For example, during the separation of minerals from their ores in the mineral processing industry, conservation of total mineral mass along with accurate prediction of the mineral fragments are truly essential. Similarly in different industries involving communition of particles, conservation of hypervolume and the zeroth moment gets basic priority. Occasionally, it is observed that certain choices of fragmentation kinetics lead to the breakdown of the mass conservative property of the particulate system [3,32]. The phenomena leading to the mass-loss from the system during fragmentation process is well known in the literature as 'shattering'. Thus, in the above mentioned industrial sectors shattering phenomenon is not at all accepted. On the other hand, a robust numerical model which estimates the physically relevant moments with good accuracy is highly acknowledged. Here, robustness of a scheme is determined by its (a) ability (straightforward or not) to get multidimensional extension, and (b) applicability on different nonuniform daughter distribution functions.
During a multidimensional fragmentation event, the internal physics of the breakage function plays a significant role in determining moments which should be conserved. This fact is well explained by [4], using schematic diagrams. Here, let us consider two example problems consisting of two properties to understand the concept.
(i) Firstly, we consider a binary breakup event of a fragment (y 1 , y 2 ) into (y 1 − x 1 , y 2 − x 2 ) and (x 1 , x 2 ). The properties x 1 and x 2 are chosen from the intervals [0, y 1 ] and [0, y 2 ], respectively. We assume that the variables y 1 and y 2 represent mass and energy, respectively. Let the kinetic kernels be given by S(y 1 , y 2 ) = 1, and b(x 1 , x 2 |y 1 , y 2 ) = 2 y 1 y 2 . (5) In this system, b is a binary breakage function and it satisfies Thus both the properties are conserved at each breakup event, and the moment equation for this problem is written as which can be easily solved to get the exact moments Thus relation (6) indicates that only the first-order moments are conserved.
In this case, the variables y 1 and y 2 represent particle length and width, respectively along two rectangular directions. For this problem, b is a multiple breakage function producing four particle fragments each of which undergoes further fragmentation equally likely. It is observed that b satisfies y1 0 y2 0 x 1 x 2 b(x 1 , x 2 |y 1 , y 2 ) dx 2 dx 1 = y 1 y 2 , and the exact moments are obtained as Unlike the previous example, the area of the particles represented by the first-cross moment M 1,1 (t) is conserved here, but it fails to conserve the first-order moments. The above mentioned examples support the fact that conservation of the moments depend upon the choice of the breakage function b( x| y). Therefore defining φ( x) := d r=1 x r , the breakage function b( x| y) requires to satisfy the relation in order to obey the first-order moment conservation law, during per fragmentation event in d−dimensions.
Similarly by definingφ( x) := d r=1 x r , the hypervolume conservation law is obeyed when b follows the relation y 0φ ( x)b( x| y) d y =φ( y).
The advances in high-speed computing have attracted many researchers to compute and simulate different particulate events. Let us have a brief review of the popular numerical methods developed in recent years to approximate the multidimensional PBEs. Primarily in most of the articles, authors have adopted different methodologies to design numerical schemes approximating the aggregation problems. In the literature, the sectional methods by [17,18] (cell average technique), [20,21,36] (fixed pivot technique), the method of higher-order moment-conserving classes by [5,6], method of moments by [24,27,39], by Monte-Carlo simulations [15,27,28,34], finite volume methods by [11,23] are well recognized because of their efficiency to predict different moments with good accuracy. However, the consideration of multidimensional fragmentation is limited to the works of [5,6,21], where the authors have designed their respective schemes to approximate the two-dimensional coupled aggregation-fragmentation equations. The work of [21] solves the fragmentation problems for uniform daughter distribution function. It is very difficult to design their schemes for daughter distribution functions which are nonuniform in nature. Moreover, extension of the [21] scheme for three or higher dimensions can be treated as a new research problem. On the other hand, the works [5] and [6] are quite similar. The method was initially proposed by [5] and then applied to solve a physical problem in [6]. In this case also, a unified formulation of [5] to solve generalized multidimensional fragmentation problems is difficult, and can be treated as the scope of future research. In this regard, the schemes based on finite volume methods are more adaptable for solving multidimensional problems. Unlike the above methods, a thorough reformulation of the finite volume scheme is not required to extend it in higher dimensions [11,13,31]. To our knowledge, approximation of multidimensional breakage problems using finite volume methods are not available in the literature till date. In this regard, an efficient numerical scheme approximating the multidimensional fragmentation models is in high demand in several industries. Here, efficiency of a numerical model is assessed upon its robustness and ability to predict different physical properties of the particulate system. Therefore we propose two finite volume schemes, designed to solve the generalized multidimensional fragmentation problems over a rectangular discretized domain. The one dimensional form of the scheme is introduced by [29]. The new schemes are formulated to predict the physically significant moments namely, the zeroth moment, the first-order moments and the d−th order moment representing particle hypervolume with good accuracy. The key feature of the schemes is that suitably defined weight functions control the number of conserved moments. Our objective is two-fold, (i) conserve the total particle properties represented by different first moments in a general multidimensional system, and (ii) conserve particle hypervolume, when each particle component represents size. Additionally in both the cases, we need to get an efficient estimation of the total particle number present in the system.
In the following section, we present a detailed discussion on the development of the proposed models. Some important follow through observations are also discussed in Section 2. Later in Section 3, we validate the efficiency of the proposed models by considering sveral test problems with two and three particle properties.
2. The multidimensional schemes. In order to perform numerical computations, we set a finite range of the computational domain. Let X := (X 1 , X 2 , . . . , X d ), where each of X r are sufficiently large positive real numbers, and hence define Ξ := d r=1 ]0, X r ] to be the truncated d−dimensional Cartesian product space. Under this consideration, the truncated equation (1) is written as . . , I with the general convention x 1r−1/2 = 0 and x Ir+1/2 = X r for all r = 1, 2, . . . , d. Thus, one can relate Ξ i to an interval in one-dimensional space, a rectangular area in two-dimensions, a cuboidal space in three-dimensions, so on and so forth. Also, let V i be the volume of the cell Ξ i , and the average particle density in each of the subspace Ξ i is defined by Let us consider , for all r = 1, 2, . . . , d and i = 1, 2, . . . , I. Letf i denotes the numerical approximation of f i over each Ξ i and S i = S( x i ). Therefore integrating (11) over each Ξ i , we obtain the following semi-discrete formulation where For brevity, we call the formulation (13) as Scheme −0. Let us now validate Scheme−0 over the example (i), from the previous section, with two internal parameters. Therefore, setting S(y 1 , y 2 ) = 1, b(x 1 , x 2 |y 1 , y 2 ) = 2 y1y2 , we compute zeroth and first moments as predicted by Scheme−0 using ODE45 solver in Matlab over a 15 × 15 rectangular mesh. Figure 1 shows that the prediction of the normalized zeroth moment by Scheme−0 is in good agreement with the exact results. Here, normalization of moments is done by dividing the values of the moments at different times by the initial value of the moment. However, Scheme−0 fails to conserve the total content of each component. In this regard, we calculate the discrete zeroth and the first moment obtained from the formulation (13). Taking sum over i on both sides of (13), we get Changing the order of the sums and simplifying, we get It can be observed that relation (15) is basically the discrete analogy of the time derivative of the continuous zeroth moment (4). Thus, a scheme which obeys (15) is expected to predict the zeroth moment with good accuracy. Let us now calculate the total discrete first moment obtained from Scheme−0. In this regard, let us define the function Φ( x i ) as the sum of all particle size representatives x ir , that is, Φ( x i ) := d r=1 x ir . Therefore, multiplying both sides by Φ( x i ) and taking sum over i, we get Proceeding in a similar fashion as done above, we can write In this case, b( x| y) satisfies the first-order moment conservation law (9). Therefore, we get Hence, It is to be noted that a scheme follows the discrete first-order moment conservation law whenever However this criterion is not true for the Scheme−0, except the trivial situation where either S or b or both are zero. Therefore in the following study, our intention to rectify Scheme−0 such that the first-order moment conservation law is obeyed.
2.1. Conservation of first-order moments. Let us first consider general multidimensional fragmentation event, where the components of x represent different particle properties like mass, enthalpy etc. Therefore to conserve the total firstorder moments, the Scheme−0 is redefined by multiplying a weight function Θ i with the second term in the right hand side, as follows with Since the formulation (17) contains one weight function, so we name it Scheme−1a for our future reference. Let us now state and prove the following proposition.
Proposition 1. The discrete system (17) is in agrement with relation (16), that is, Scheme−1a conserves the total first-order moments of the system.
Proof. The proposition can easily be proved by performing similar calculations as done in the preceding section.
The above proposition suggests that Scheme−1a obeys the first-order moment conserving criterion of the multidimensional system. However, Scheme−1a does not follows the discrete formulation (15). Taking sum over i on both sides of (17), we get Therefore, from the previous observation we can say that Scheme−1a cannot predict the total particle number with good accuracy. Since, our aim is to obtain a scheme which conserves the total first-order moments as well as, predicts zeroth moment efficiently. Therefore, the formulation (13) is further modified by introducing two weights in the right hand side, to get where, Considering that formulation (19) involves two weight functions, we call it Scheme −2a. Basically the weight functions, Ψ b k and Ψ d i control particle distribution in order to conserve the zeroth and first moments. Let us now state and prove the following proposition.
Proposition 2. The numerical model (19) conserves the total first-order moments (16), and also obeys the discrete zeroth number prediction given by the relation (15), whenever the weight functions Ψ b k and Ψ d k satisfy the formulations of (20). Proof. We compute the total temporal change of the discrete first moments, and get .
The temporal change of discrete zeroth moment is computed as follows, Thus Scheme−2a (19) conserves the first-order moments, and also follows the temporal change of discrete zeroth moment (15).

2.2.
Conservation of hypervolume. We now consider a fragmentation event, where each component of x represents particle size. As mentioned in Section 1, the hypervolume of the particles plays an important role in several real life occasions, as it represents the geometry of the particles. Therefore depending upon the situation, the hypervolume also needs to get conserved, that is, a numerical model should obey where,Φ In this part of our study, we redefine the weight functions such that the new set of schemes conserve particle hypervolume. Moreover, during hypervolume conservation per fragmentation event, the breakage function b( x| y) should satisfy (10).
Therefore, Scheme−1b is written as with the modified weightΘ and Scheme−2b is written as having the modified weights It can be easily proved that both the Scheme−1b (22) and Scheme−2b (24) follows hypervolume conservation law (21). Additionally, Scheme−2b (24) follows the discrete number formulation given by (15). The proofs of these claims bears quite similarity with the previous ones, therefore we omit them here. Therefore, in this section we systematically propose two sets of finite volume schemes. The first set consisting of Scheme−1a [eqs (17), (18)] and Scheme−2a [eqs (19), (20)] are basically designed to conserve the sum of the first-order moments (16). Theoretically, Scheme−2a estimates the evolution of total number of particles with high accuracy. In the second set, two models Scheme−1b [eqs (22), (23)] and Scheme−2b [eqs (24), (25)] are designed to conserve the particle hypervolume (21) during a breakage event. In this case also, the Scheme−2b is expected to estimate the zeroth moment with good accuracy. Based on their ability to predict two physically significant moments, Scheme−2a and Scheme−2b are the desired models which should be used to solve multidimensional fragmentation events. Remark 1. It is to be noted that the weight functions corresponding to the firstorder moment and hypervolume conservation are interrelated. Basically, it depends upon the choice of the breakage function b( x| y). Therefore, the factor Φ ( x i ) representing the sum of the pivots in the first-order moment conservative models is simply replaced by the factorΦ ( x i ) representing the product of the pivots in the hypervolume conservative models.
3. Test cases and numerical details. In this section, we validate the efficiency of Scheme−1a (17), Scheme−2a (19) and Scheme−1b (22), Scheme−2b (24) by considering eight test problems defined by two particle properties, and two test problems with three particle properties. As mentioned in Section 1, three dimensional extension of the schemes of [21] and [5] are very difficult task and can be treated as a scope of future research. Therefore, to maintain uniformity throughout the two and three dimensional test problems presented in the current paper, we have not included the schemes of [5,21] for comparison. Moreover, computing the particle property distribution and its moments is a challenging task for multidimensional breakage equations, as the literature lacks to provide the exact solutions in closed form. So, it is not always possible to compare the numerical number density against the exact value. Instead, we consider certain sample problems with physically relevant daughter distribution function for which the moments can be calculated exactly. The example problems are gathered from the articles of [4,5,16,26].
In the following study, the two-dimensional test problems are grouped into two sets. First set consists of four test problems which are chosen to conserve the first-order moments (16). The second set with another four problems which conserve hypervolume of particle property distribution (21). Moreover, for each of the above mentioned test problems, we choose both uniform and nonuniform breakage functions to validate the efficiency of the proposed schemes.
Initially, we consider a size-independent constant selection function S 1 (x 1 , x 2 ) = 1, along with the following daughter distribution functions b 1 (x 1 , x 2 |y 1 , y 2 ) = 2 y 1 y 2 , and b 2 (x 1 , x 2 |y 1 , y 2 ) = 2δ Here, b 1 represents the uniform distribution of daughter particles, and b 2 represents a nonuniform symmetric daughter distribution function. Both the functions b 1 and b 2 satisfy relation (9). We choose a mono-dispersed initial data for test case 1. The moment function for b 1 has already been calculated in relation (6), and for b 2 it is calculated in [5] and is written as For mono-dispersed initial data, we get M k,l (0) = 1.
In the second instance, we consider a size dependent selection function S 2 (x 1 , x 2 ) = x 1 + x 2 . Both the breakage functions b 1 and b 2 mentioned above, are recalled here along with the mono-dispersed initial data. In this case, only the moments M 0,0 (t), M 1,0 (t) and M 0,1 (t) can be calculated exactly. Interestingly, the moment functions obtained for both the daughter distribution function are same, and are given by M 1,0 (t) = M 0,1 (t) = 1 and M 0,0 (t) = 1 + 2t.
Next to it, four examples are considered where the daughter distribution functions are designed to conserve particle hypervolume, that is, the breakage functions satisfy relation (10). In this regard, we first consider a size-independent selection function S 3 (x 1 , x 2 ) = 1. The following breakage functions are chosen for numerical evaluation Physical interpretation of the breakage function b 3 can be found in the article of [16], and that of b 4 in [26]. Relation (8) gives the closed form of exact moments with kernels S 3 and b 3 . Similarly, for S 3 and b 4 the following exact moments are calculated as Similarly, the moments for breakage function b 4 and selection function S 4 are calculated as M 1,1 (t) = 1 and M 0,0 (t) = 1 + t.
In the later part of this study, we consider one example from each of the above mentioned problem sets in their three-dimensional variant. A three-dimensional mono-dispersed initial data is also considered to support the problems. Here, two test problems are defined with the following kernels (a) y1y2y3 . The exact moments corresponding to the the kernels S 5 , b 5 are calculated as M 1,0,0 (t) = M 1,0,1 (t) = M 0,0,1 (t) = 1 and M 0,0,0 (t) = 1 + 3t.
In the following tables, we summarize all the test problems mentioned above for a smooth reading of the subsequent section. From Table 1, we find that a general closed form of the higher order exact moments can only be found for the test case 1, 2 and 5. For rest of the sample problems only first few moments can be calculated exactly. It can be observed that test case 9 is a straightforward extension of test case 4 with symmetric breakage function, and the test case 10 is an extension of test case 7 with uniform breakage function in three-dimensions. Table 2. Summary of the selected test problems in three dimensions.
The efficiency of the proposed schemes is analyzed both qualitatively and quantitatively. The qualitative comparison includes the graphical representation of the moment functions, as presented in the literature by [13,31]. To obtain a systematic plot, the moments are sorted in decreasing order of their exact values. On the other hand, weighted relative errors are evaluated at different times to determine the accuracy of the moment functions quantitatively. A general measure of the weighted relative error for two-dimensional problems is given as [31] .
Furthermore to obtain a clear visibility of different markers, all the figures are plotted in linear scale along horizontal-axis and in logarithmic scales along the vertical-axis. All the computations in this part are carried till T = 5 in a standard computer with i3 (2.2 GHz r.p.m.) processor and 3GB RAM. We use Ode45 solver in Matlab to compute the system of differential equations. The CPU usage time required for solving the test problems are also included. 4. Results and discussion.

Conservation of first-order moment. Size independent selection function:
In Figure 2, graphical validation of the dimensionless moments M 0,0 (t), M 1,0 (t) and M 1,1 (t) are presented for the test case 1 (Figure 2a) and test case 2 ( Figure  2b). The quantitative analysis of the proposed schemes is presented in the tables 3, 4, 5 and 6. The tables 3 and 4 contain weighted relative errors accumulated to estimate the three moments M 0,0 (t), M 1,0 (t) and M 1,1 (t) using Scheme−1a (17) and Scheme−2a (19). Table 3 represents the errors for the test case 1, and Table 4 corresponds the same for the test case 2. In some occasions, higher order moments play important roles. Therefore for a complete quantitative analysis, we present the weighted relative errors accumulated by the schemes to estimate higher order moments over different grids at the end time t = 5. Similar, presentation of the errors can be found in the article of [5]. The data for test case 1 are given in Table  5 and that of test case 2 in Table 6. From Figure 2, it is observed that Scheme−2a conserves the first moment M 1,0 (t) and predicts the zeroth moment M 0,0 (t) with high accuracy. Thus rectifying the drawbacks of Scheme−0 (13). It also produces a good estimation the first cross moment M 1,1 (t). The estimation of the moment M 1,1 (t) improves with further refinement of the meshes. On the other hand, Scheme−1a conserves the first moment, but the prediction of the zeroth moment M 0,0 (t) and the cross moment M 1,1 (t) is quite poor. This observation is also validated quantitatively by the weighted error Table 3 and Table 4.  Table 3. Relative error for the weighted moments at different times for the test case 1.
The Scheme−1a conserves the sum of all first-order moments. The other moments, namely M 0,0 (t) and M 1,1 (t) are not conserved by Scheme−1a. From the data given in Table 3 and Table 4, we observe that the error accumulated to estimate the conserved first-order moments is negligible as compared to the errors accumulated to estimate the non-conserved moments M 0,0 (t) and M 1,1 (t). Similarly, Scheme−2a efficiently estimates the zeroth and first-order moments, but fails to predict M 1,1 (t) with good accuracy.  To solve test case 1 and test case 2, the CPU times taken by Scheme−1a is approximately 1s, and that by Scheme−2a is 2s (nearly).
From Table 5 and Table 6, we observe that the estimation of the higher order moments are not as good as the zeroth and the first order moments despite of several refinement of the meshes.
Size dependent selection function: Here, we discuss the efficiency of the proposed schemes for solving test case 3 and test case 4. Numerical values of the moments M 0,0 (t) and M 1,0 (t) are plotted against their exact values in Figure 3 and the weighted errors are calculated in Table 7 and Table 8. Figure 3, Table 7 and Table  8 indicate that Scheme−2a is highly accurate in estimating the moments compared to Scheme−1a. The exact higher order moments can not be obtained for this problem.

4.2.
Conservation of hypervolume. Size independent selection function: In Figure 4, we plot the numerical values of the moments M 0,0 (t), M 1,0 (t) and M 1,1 (t) obtained from Scheme−1b (22) and Scheme−2b (24) against their exact values. On the other hand, Table 10 and Table 11 Table 6. Relative error for higher order weighted moments using different computational grids for the test case 2 at t = 5.       obtained from test case 5 and test case 6, respectively. However, the higher order moments can be calculated exactly only for the test case 5. Therefore, the quantitative error analysis for Scheme−1b and Scheme−2b for higher moments are given in Table 12. As expected, Scheme−2b produces highly improved estimation of the moments M 0,0 (t) and M 1,1 (t) over Scheme−1b. From Figure 4a, it is observed that Scheme−2b overpredicts the first moment M 1,0 (t) for the kernel b 3 . However, the accuracy increases considerably for test problem 6 over the same number of meshes [ Figure 4b].
In the above tables, we observe that Scheme−1b conserves the moment M 1,1 (t), whereas the other non-conserved moments are poorly predicted. Similarly, Scheme− 2b accumulates nearly negligible error to predict the conserved M 0,0 (t) and M 1,1 (t) moments. However, it estimates the other non-conserved moments poorly. The CPU time taken by Scheme−1b to solve test cases 5 and 6 is approximately 1s each. Similarly, time taken by Scheme−2b is nearly 2s for both the problems.
Size dependent selection function: The Figure 5, Table 13 and Table 14 validate that the Scheme−2b predicts the zeroth and first cross moments of the particles property distribution with high accuracy. Moreover, both the schemes take less than 1s to perform the computations.
The CPU usage time (in seconds) taken to solve test cases 7 and 8 is nearly 1s.

4.3.
Test cases in three-dimensions. In this part, Figure 6a represents the qualitative efficiency of the proposed schemes to predict different moments. Conservation of the total first-order moments is illustrated in Figure 6a. Similarly, the conservation of the third-order moment corresponding to particle hypervolume is shown in Figure 6b. Additionally, Table 15 and Table 16 represent the weighted relative errors accumulated by the two schemes at different times for the test cases 9 and 10, respectively. In both the cases, we find that the schemes with two weight functions outperforms the schemes with one weight function to estimate the zeroth moment representing the total particle number. Interestingly, Table 17 depicts that both the schemes consume low CPU time to perform the computations for a three dimensional model. Other higher order moments cannot be evaluated exactly.
Here, Scheme−1a is the three-dimensional extension of the models (17), (18

5.
Conclusions. In this article, we have introduced two sets of finite volume schemes approximating the multidimensional fragmentation problems. The schemes are designed to estimate physically important moment functions of the particle property distribution with high accuracy. It is observed that the efficiency of the schemes with two weight functions in predicting the zeroth moment is considerably higher compared to the single weighted scheme. Therefore as a natural selection, Scheme−2a (19) and Scheme−2b (24) are our preferred models to approximate the solutions of multidimensional fragmentation PBEs. However, both the schemes have simple mathematical formulations and are robust to compute on nonuniform meshes. Moreover, the above schemes consumes a very low CPU usage time to solve multidimensional fragmentation problems. Thus one can easily compute the proposed models in a standard PC and solve fragmentation problems in any dimension/s. The efficiency of the schemes are validated over several test problems.