*Article* **Conservative Finite Volume Schemes for Multidimensional Fragmentation Problems**

**Jitraj Saha <sup>1</sup> and Andreas Bück 2,\***


**Abstract:** In this article, a new numerical scheme for the solution of the multidimensional fragmentation problem is presented. It is the first that uses the conservative form of the multidimensional problem. The idea to apply the finite volume scheme for solving one-dimensional linear fragmentation problems is extended over a generalized multidimensional setup. The derivation is given in detail for two-dimensional and three-dimensional problems; an outline for the extension to higher dimensions is also presented. Additionally, the existing one-dimensional finite volume scheme for solving conservative one-dimensional multi-fragmentation equation is extended to solve multidimensional problems. The accuracy and efficiency of both proposed schemes is analyzed for several test problems.

**Keywords:** conservative formulation; multidimensional fragmentation equation; weight functions; finite volume scheme

**Citation:** Saha, J.; Bück, A. Finite Volumes Schemes for Multidimensional Fragmentation Problems. *Mathematics* **2021**, *9*, 635. https://doi.org/10.3390/math9060635

Academic Editor: Lucas Jódar

Received: 8 February 2021 Accepted: 9 March 2021 Published: 17 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

#### **1. Introduction**

Fragmentation, breakage or attrition, describe processes in which a single object is separated into at least two new objects. The reasons for breakage can be manifold but are often linked to some kind of stress exerted on the object, for instance thermal stress from heating and rapid cooling (or vice versa and cyclically)—a natural process specifically observed in deserts, leading to disintegration of rocks; or mechanical stress, applied for millenia in the process of grain milling. Nowadays, fragmentation plays a key role in several industrial sectors like mineral processing (e.g., comminution of ores [1–4]), reaction engineering (e.g., break-up of bubbles in reacting bubble columns for separation processes [5–8] or steel-casting [9]) or pharmaceutical industries (e.g., milling of active pharmaceutical ingredients to increase their solubility and uptake capacity in the human or animal body [10–12]).

Many objects, e.g., particles, bubbles or even rain drops, consist of different components resulting in their anisotropic structure. The probability of fragmentation upon stress therefore depends on the distribution of the components within the objects, i.e., each component adds an independent dimension to the fragmentation problem. In an attempt to describe these complex processes and make them accessible for model- and knowledge-based process design, optimization and control, multidimensional fragmentation equations have been proposed and used in different fields of application, see, for instance, the works [13–17].

Theoretical aspects on the existence of scaling solutions and their behavior at the onset of "shattering" transition have been discussed for instance in the works of [18–21]. Fragmentation models are particularly challenging as they consist of partial-integro differential equations as will be shown in the following. Analytical results are scarce and often of very limited practical relevance, strongly motivating the development of numerical methods for approximation of the solution to (multidimensional) fragmentation problems.

As a prototype, consider the conservative formulation of the multiple fragmentation equation given by [22,23]: The initial value problem for *t* ≥ 0 is formulated as

$$\frac{\partial \mathbf{g}(t, \mathbf{x})}{\partial t} = \frac{\partial \mathcal{H}(t, \mathbf{x})}{\partial \mathbf{x}}, \quad \text{where} \quad \mathbf{x} \in \mathbb{R}^+ := (0, \infty), \tag{1}$$

with the initial condition

$$\log(\mathbf{x},0) = \mathfrak{g}\_0(\mathbf{x}) \quad (\geq 0), \quad \text{for} \quad \mathbf{x} \in \mathbb{R}^+.\tag{2}$$

The *flux function* <sup>H</sup>(*t*, *<sup>x</sup>*) is defined by

$$\mathcal{H}(t, \mathbf{x}) := \int\_{x}^{\infty} \int\_{0}^{\mathbf{x}} \frac{\mathbf{u}}{\upsilon} b(\mathbf{u}, \upsilon) S(\upsilon) g(t, \upsilon) \mathbf{d}u d\upsilon, \quad \mathbf{x} \in \mathbb{R}^{+}. \tag{3}$$

In Equation (1), the internal variables *x* and *t* denote the particle property and the time component, respectively. On the left hand side, the function *<sup>g</sup>*(*t*, *<sup>x</sup>*) is defined by *<sup>g</sup>*(*t*, *<sup>x</sup>*) := *x f*(*t*, *<sup>x</sup>*), where *<sup>f</sup>*(*t*, *<sup>x</sup>*) denotes the distribution of particle volume *<sup>x</sup>* in a system at time *t*. The rate of selection of an *x*-volume cluster to undergo breakage is denoted by *<sup>S</sup>*(*x*), and the distribution of daughter particles *<sup>y</sup>* due to the breakage of large particle *<sup>x</sup>* is denoted by *<sup>b</sup>*(*y*, *<sup>x</sup>*). The breakage function *<sup>b</sup>*(*y*, *<sup>x</sup>*) satisfies the following relations:

$$\int\_0^\infty b(y, \mathbf{x}) \mathbf{d}y = \nu(\mathbf{x}), \quad \text{and} \quad \int\_0^\mathbf{x} y b(y, \mathbf{x}) \mathbf{d}y = \mathbf{x}. \tag{4}$$

The first relation defines that *<sup>ν</sup>*(*x*) number of fragments are produced during the breakup of a large *x*-cluster, and the second relation defines that the total volume of the daughter *y*-clusters is exactly the same as the volume of the mother *x*-cluster. Note that the formulation (1) is well-known in the literature as the volume conservative form. Integration of Equation (1) over the volume variable *x* from 0 to ∞, with the help of relation (4), yields

$$\frac{d}{dt} \int\_0^\infty g(t, x) dx = 0. \tag{5}$$

It should be noted that Equation (1) is a first-order hyperbolic, initial value partial differential equation. In this regard, the representation (1) gathered importance because the divergent nature allows the model to obey the volume conservation laws. The coefficient *S* belongs to *L*<sup>∞</sup> *loc*([0, <sup>∞</sup>)) and *<sup>g</sup>*0, *<sup>b</sup>* <sup>∈</sup> *<sup>L</sup>*1((0, <sup>∞</sup>)) <sup>∩</sup> *<sup>L</sup>*1((0, <sup>∞</sup>), *<sup>x</sup>*d*x*). Here and below, the notation *<sup>L</sup>*1(R+, *<sup>x</sup>*d*x*) stands for the space of the Lebesgue measurable real-valued functions on R<sup>+</sup> which are integrable with respect to the measure *<sup>x</sup>*d*x*.

In most of the previous studies it is assumed that a single parameter, which is usually volume, mass or size of the particle, is sufficient to describe the particle property (readers can refer to [24] for further details). However, a single parameter is not always sufficient to describe various physical systems. For example, fragment mass distribution obtained by crushing gypsum or glass depends on the initial geometry of the particles. On the other hand, the degradation of polyelectrolyte may depend upon both their mass and excitation (or kinetic) energy. Therefore, the fragmentation dynamics need to be represented by including additional variables to the mathematical model. These variables are equivalently classified as the degrees of freedom of the dynamical system and hence, the multidimensional formulation of the fragmentation equations becomes necessary to represent such cases. The purpose of this article is to take in account more than one particle property and present an efficient numerical model which estimates them with high accuracy. In particular, we present the mathematical representations of two-dimensional and threedimensional volume conservative linear fragmentation equations. Further extension of the mathematical formulation can be done in a similar manner.

For the population balance models, the moment functions of the particle property distribution play a major role as some of them describe a significant physical property of the system. Therefore, before we proceed further, let us first gather some important information about the moment functions in a generalized multidimensional setup.

#### *1.1. Moment Functions*

Let **<sup>x</sup>** := {*x*1, *<sup>x</sup>*2, ... , *xn*}, with *xi*−s representing different particle properties like, mass, entropy, moisture content, shape factor, etc. and thus the function *<sup>f</sup>*(*t*, **<sup>x</sup>**) denotes the distribution of particle property **x** at some instance *t*. The formal definition of the moment functions for a general *n*-dimensional population balance problem is written as follows:

$$\mathcal{M}\_{p\_1, p\_2, \dots, p\_n}(t) := \int\_0^\infty \left( \prod\_{r=1}^n \mathbf{x}\_r^{p\_r} \right) f(t, \mathbf{x}) d\mathbf{x}\_r \tag{6}$$

where the integrations are defined as

$$\int\_0^\infty \quad (\cdot) \mathbf{dx} := \underbrace{\int\_0^\infty \mathbf{dx}\_1 \cdots \int\_0^\infty}\_{n-\text{times}} \quad (\cdot) \mathbf{dx}\_n \cdots$$

In Equation (6), *p*1, *p*2, ... , *pn* are nonnegative integers. As mentioned earlier, the moment functions play an important role to define various physical properties of the system. Like the zeroth moment, <sup>M</sup>0,0,...,0(*t*) defines the total number of particles present in the system. The first-order moment <sup>M</sup>0,...,1,...,0(*t*) (1 is the *<sup>r</sup>*th position) denotes the total content of the *xr*th component in the system, which can be equivalently represented as the total volume of *xr*th property. Hence, for a multidimensional system, the volume conservation of the system can be defined as the total conservation of all the first-order moments taken together. Thus, defining *<sup>φ</sup>*(**x**) := *n* ∑ *xr*, the volume conservation law for

*r*=1

the *n*-dimensional system is expressed as

$$\frac{d}{dt} \int\_0^\infty \phi(\mathbf{x}) f(t, \mathbf{x}) d\mathbf{x} = 0. \tag{7}$$

Furthermore, the *<sup>n</sup>*-th order cross moment is defined by <sup>M</sup>1,...,1,...,1(*t*) and it represents the particle geometry or hypervolume. Therefore, to preserve the initial geometry of the particles, we need to preserve the cross moments; hence, the hypervolume preservation law is written as

$$\frac{d}{dt}\int\_{0}^{\infty}\psi(\mathbf{x})f(t,\mathbf{x})d\mathbf{x}=0,\tag{8}$$

where *<sup>ψ</sup>*(**x**) :<sup>=</sup> *n* ∏*r*=1 *xr*.

Similarly, other higher order moments can be defined using the formulation (6), and depending upon the problem they may correlate to some physical properties of the system. For example, in a pipeline flow for the transport of natural gas from seabed, the breakage of hydrate particles often takes place. In this event, if the first moment <sup>M</sup>1,...,0(*t*) is proportional to the mean radius of the hydrate particle, then the corresponding second order moment <sup>M</sup>2,...,0(*t*) and third order moment <sup>M</sup>3,...,0(*t*) are proportional to the total area and the volume concentration of the hydrate particles, respectively. In general, only the zeroth, first-order and the cross-moments bear the same meaning for any population balance models.However, it should not be misunderstood that higher order moments should always correspond to certain physical characteristics.

In the literature, a limited number of articles are dedicated to the numerical study of multidimensional fragmentation events, and therefore several aspects of study still remain unexplored. The articles of [25–29] discuss the development of different numerical schemes to approximate the fragmentation problems. To note that, unlike the methods, e.g., cell average technique, fixed pivot techniques, method of moments, etc., the finite volume methods have gained popularity because the latter are robust to be applied on a multidimensional setup. Moreover, the underlying stencil of the finite volume scheme is simpler, and easy to compute (the readers can refer to the articles of [23,29] for further details on the computational advantage of finite volume schemes).

The article is organized in the following manner. In the next section, we present the mathematical representations of the continuous two- and three-dimensional equations. In this regard, the three-dimensional model is represented using vector notation, which will also provide an outline to extend the equations into further higher dimensions. In Section 3, step-by-step derivation of the numerical schemes are presented. An interesting outcome of this presentation includes the multidimensional extension of the finite volume scheme presented in [23]. Section 4 contains the numerical validation of the proposed models over some standard empirical test problems. Finally some conclusions and a summary of the work are presented.

#### **2. Continuous Equations in Two- and Three-Dimensions**

#### *2.1. Conservative Formulations in Two-Dimensions*

In Equation (1), the variable *x* represents a single particle property which can be considered as the particle volume. Therefore, the first moment always corresponds to the total volume of the particle in the system, and hence Equation (1) is simply coined as the volume-conservative model. However, the representation is not that simple in the case of a multidimensional fragmentation event. Depending on the definition of volume and hypervolume, the mathematical model changes, and thus we get two different mathematical equations representing the two conservative formulations in the multidimensional setup. For example, consider two independent particle properties kinetic energy and moisture content that are defined by the variables *<sup>x</sup>* and *<sup>y</sup>*, respectively and we set **<sup>x</sup>** := (*x*, *<sup>y</sup>*). Then *<sup>f</sup>*(*t*, **<sup>x</sup>**) is the two-dimensional particle properties distribution function at time *<sup>t</sup>*. Now referring to the Equations (7) and (8), the solutions corresponding to the volume-conservative and hypervolume conservative formulations are defined as *<sup>n</sup>*(*t*, **<sup>x</sup>**) := *<sup>φ</sup>*(**x**)*f*(*t*, **<sup>x</sup>**) and *<sup>m</sup>*(*t*, **<sup>x</sup>**) := *<sup>ψ</sup>*(**x**)*f*(*t*, **<sup>x</sup>**), respectively.

In accordance with the above definition, the two-dimensional or bivariate volumeconservative fragmentation equation is written as,

$$\frac{\partial \mathfrak{n}(t, \mathbf{x})}{\partial t} = \frac{\partial \mathcal{F}(t, \mathbf{x})}{\partial \mathbf{x}} + \frac{\partial \mathcal{G}(t, \mathbf{x})}{\partial y} - \frac{\partial^2 \mathcal{H}(t, \mathbf{x})}{\partial \mathbf{x} \partial y},\tag{9}$$

with the initial data

$$m(0, \mathbf{x}) = m\_0(\mathbf{x}) \ge 0, \quad \text{for all} \quad \mathbf{x} > 0. \tag{10}$$

Here, the functions F, G and H denote the flux flow at the cell boundaries. In this regard, we first define the following notations to be used for defining the fluxes. Let **<sup>u</sup>** := (*u*, *<sup>v</sup>*), := (, *<sup>ξ</sup>*), then

$$\int\_{\mathbf{x}}^{\infty} (\cdot) \mathbf{d} \mathbf{u} := \int\_{\mathbf{x}}^{\infty} \int\_{y}^{\infty} (\cdot) \mathbf{d} \mathbf{v} \mathbf{d} \mu, \quad \text{and} \quad \int\_{0}^{\mathbf{x}} (\cdot) \mathbf{d} \mathbf{u} := \int\_{0}^{\mathbf{x}} \int\_{0}^{y} (\cdot) \mathbf{d} \mathbf{v} \mathbf{d} \mu.$$

With the help of the above notations, the flux functions are defined as follows:

$$\mathcal{F}(t, \mathbf{x}) := \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}} \frac{\phi(\boldsymbol{\varepsilon}, \boldsymbol{y})}{\phi(\mathbf{u})} b(\boldsymbol{\varepsilon}, \boldsymbol{y} | \mathbf{u}) S(\mathbf{u}) n(t, \mathbf{u}) d\boldsymbol{\varepsilon} d\mathbf{u},\tag{11}$$

$$\mathcal{G}(t, \mathbf{x}) := \int\_{\mathbf{x}}^{\infty} \int\_{0}^{y} \frac{\phi(\mathbf{x}, \frac{\mathbf{x}}{\sigma})}{\phi(\mathbf{u})} b(\mathbf{x}, \mathbf{\tilde{g}} | \mathbf{u}) S(\mathbf{u}) n(t, \mathbf{u}) d\mathbf{\tilde{g}} d\mathbf{u},\tag{12}$$

and

$$\mathcal{H}(t, \mathbf{x}) := \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}} \frac{\phi(\mathbf{e})}{\phi(\mathbf{u})} b(\mathbf{e}|\mathbf{u}) S(\mathbf{u}) n(t, \mathbf{u}) \mathrm{d}\mathbf{e} d\mathbf{u}.\tag{13}$$

In the above expressions, *<sup>S</sup>*(**x**) is the selection function which defines the rate at which particles of properties **<sup>x</sup>** to undergo further fragmentation, and the breakage function *<sup>b</sup>*(|**u**) corresponds to the distribution of daughter fragments formed due to the fragmentation of **u**-cluster. In the multidimensional fragmentation setup, the breakage function *b* plays a key role to govern the system to obey either volume-conservation (7) or hypervolume conservation (8) laws. For the volume conservative formulation, it is assumed that the breakage function *<sup>b</sup>*(|**u**) should satisfy

$$\int\_0^\mathbf{u} \phi(\boldsymbol{\varepsilon}) b(\boldsymbol{\varepsilon}|\mathbf{u}) d\boldsymbol{\varepsilon} = \phi(\mathbf{u}).\tag{14}$$

The relation (14) is significant as it controls the system to follow volume conservation property (7). Therefore, with the above assumption it can easily be calculated that

$$\frac{\mathrm{d}}{\mathrm{d}t} \int\_{0}^{\infty} n(t, \mathbf{x}) \mathrm{d}\mathbf{x} = \frac{\mathrm{d}}{\mathrm{d}t} [\mathcal{M}\_{1,0}(t) + \mathcal{M}\_{0,1}(t)] = 0,\tag{15}$$

that is, the volume conservation laws are perfectly obeyed.

Note that the flux function H (13) in the bivariate Equation (9) is the straightforward extension of the flux in univariate model (1). Additionally, the bivariate model (9) contains two flux functions F (11) and G (12) as compared to its one-dimensional counterpart Equation (1). Here, F defines the distribution of the daughter particles along the *x*component, while keeping the *y*-component fixed. Similarly, the flux G is defined along *y*-component.

We now present the continuous hypervolume conservative formulation of the pure bivariate fragmentation model. It is expressed in a manner similar to the volume conservative model (9), and reads as

$$\frac{\partial m(t, \mathbf{x})}{\partial t} = \frac{\partial \vec{\mathcal{F}}(t, \mathbf{x})}{\partial x} + \frac{\partial \vec{\mathcal{G}}(t, \mathbf{x})}{\partial y} - \frac{\partial^2 \vec{\mathcal{H}}(t, \mathbf{x})}{\partial \mathbf{x} \partial y},\tag{16}$$

with the flux functions <sup>F</sup>¯, <sup>G</sup>¯ and <sup>H</sup>¯ redefined as follows

$$\mathcal{F}(t, \mathbf{x}) := \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\chi} \frac{\psi(\epsilon, y)}{\psi(\mathbf{u})} b(\epsilon, y|\mathbf{u}) S(\mathbf{u}) m(t, \mathbf{u}) d\epsilon d\mathbf{u},\tag{17}$$

$$\mathcal{G}(t, \mathbf{x}) := \int\_{\mathbf{x}}^{\infty} \int\_{0}^{y} \frac{\psi(\mathbf{x}, \underline{\xi})}{\psi(\mathbf{u})} b(\mathbf{x}, \underline{\xi} | \mathbf{u}) S(\mathbf{u}) m(t, \mathbf{u}) \, \mathrm{d}\underline{\xi} \, \mathrm{d}\mathbf{u},\tag{18}$$

and

$$\bar{\mathcal{H}}(t, \mathbf{x}) := \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\infty} \frac{\psi(\epsilon)}{\psi(\mathbf{u})} b(\epsilon | \mathbf{u}) S(\mathbf{u}) m(t, \mathbf{u}) d\epsilon d\mathbf{u}.\tag{19}$$

In this case, the breakage function satisfies condition

$$\int\_0^\mathbf{u} \psi(\boldsymbol{\epsilon}) b(\boldsymbol{\epsilon}|\mathbf{u}) d\boldsymbol{\epsilon} = \psi(\mathbf{u}),\tag{20}$$

and hence, integrating Equation (16), one can easily obtain that the hypervolume conservation law is properly obeyed, that is

$$\frac{d}{dt} \int\_0^\infty m(t, \mathbf{x}) d\mathbf{x} = \frac{\mathbf{d} \mathcal{M}\_{1,1}(t)}{\mathbf{d}t} = 0. \tag{21}$$

#### *2.2. Conservative Formulations in Three-Dimensions*

In a similar manner as discussed above, we now present the three-dimensional representations of the fragmentation equations which obey the (i) volume conservation laws, and (ii) hypervolume conservation laws. In this part, we present the mathematical model using vector notation to pave the way for higher dimensional extension.

Let the particle property distribution be written as *<sup>f</sup>*(*t*, **<sup>x</sup>**) where the vector **<sup>x</sup>** := {*x*1, *<sup>x</sup>*2, *<sup>x</sup>*3} represents different particle properties, and *<sup>n</sup>*(*t*, **<sup>x</sup>**) := *<sup>φ</sup>*(**x**)*f*(*t*, **<sup>x</sup>**). Using the extended form of all the above-mentioned notations, the three-dimensional volume conservative formulation is written as follows

$$\begin{split} \frac{\partial n(t,\mathbf{x})}{\partial t} &= \frac{\partial \mathcal{F}^{(1)}(t,\mathbf{x})}{\partial \mathbf{x}\_{1}} + \frac{\partial \mathcal{F}^{(2)}(t,\mathbf{x})}{\partial \mathbf{x}\_{2}} + \frac{\partial \mathcal{F}^{(3)}(t,\mathbf{x})}{\partial \mathbf{x}\_{3}} \\ &- \frac{\partial^{2} \mathcal{G}^{(1)}(t,\mathbf{x})}{\partial \mathbf{x}\_{2} \partial \mathbf{x}\_{3}} - \frac{\partial^{2} \mathcal{G}^{(2)}(t,\mathbf{x})}{\partial \mathbf{x}\_{1} \partial \mathbf{x}\_{3}} - \frac{\partial^{2} \mathcal{G}^{(3)}(t,\mathbf{x})}{\partial \mathbf{x}\_{1} \partial \mathbf{x}\_{2}} + \frac{\partial^{3} \mathcal{H}(t,\mathbf{x})}{\partial \mathbf{x}\_{1} \partial \mathbf{x}\_{2} \partial \mathbf{x}\_{3}}, \end{split} \tag{22}$$

with the flux flows being functions of both *t*, **x** and are defined as

$$\mathcal{F}^{(1)}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}\_1} \frac{(\mathbf{u}\_1 + \mathbf{x}\_2 + \mathbf{x}\_3)}{\phi(\mathbf{y})} b(\mathbf{u}\_1, \mathbf{x}\_2, \mathbf{x}\_3 | \mathbf{y}) S(\mathbf{y}) n(t, \mathbf{y}) d\mathbf{u}\_1 d\mathbf{y},\tag{23}$$

$$\mathcal{F}^{(2)}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}\_2} \frac{(\mathbf{x}\_1 + \mathbf{u}\_2 + \mathbf{x}\_3)}{\phi(\mathbf{y})} b(\mathbf{x}\_1, \mathbf{u}\_2, \mathbf{x}\_3 | \mathbf{y}) \mathcal{S}(\mathbf{y}) n(t, \mathbf{y}) d\mu\_2 d\mathbf{y},\tag{24}$$

$$\mathcal{F}^{(3)}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}\_{3}} \frac{(\mathbf{x}\_{1} + \mathbf{x}\_{2} + \mathbf{u}\_{3})}{\phi(\mathbf{y})} b(\mathbf{x}\_{1}, \mathbf{x}\_{2,3} \, | \, \mathbf{y}) \mathcal{S}(\mathbf{y}) n(t, \mathbf{y}) d\mathbf{u}\_{3} d\mathbf{y},\tag{25}$$

$$\mathcal{G}^{(1)}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}\_2} \int\_{0}^{\mathbf{x}\_3} \frac{(\mathbf{x}\_1 + \mathbf{u}\_2 + \mathbf{u}\_3)}{\phi(\mathbf{y})} b(\mathbf{x}\_1, \mathbf{u}\_2, \mathbf{u}\_3 | \mathbf{y}) \mathcal{S}(\mathbf{y}) n(t, \mathbf{y}) d\mathbf{u}\_2 d\mathbf{u}\_3 d\mathbf{y},\tag{26}$$

$$\mathcal{G}^{(2)}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}\_1} \int\_{0}^{\mathbf{x}\_3} \frac{(\mathbf{u}\_1 + \mathbf{x}\_2 + \mathbf{u}\_3)}{\phi(\mathbf{y})} b(\mathbf{u}\_1, \mathbf{x}\_2, \mathbf{u}\_3 | \mathbf{y}) \mathcal{S}(\mathbf{y}) n(t, \mathbf{y}) d\mathbf{u}\_3 d\mathbf{u}\_1 d\mathbf{y}\_\prime \tag{27}$$

$$\mathcal{G}^{(3)}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}\_1} \int\_{0}^{\mathbf{x}\_2} \frac{(\boldsymbol{\mu}\_1 + \boldsymbol{\mu}\_2 + \mathbf{x}\_3)}{\phi(\mathbf{y})} b(\boldsymbol{\mu}\_1, \boldsymbol{\mu}\_2, \mathbf{x}\_3 | \mathbf{y}) \mathcal{S}(\mathbf{y}) n(t, \mathbf{y}) d\boldsymbol{\mu}\_2 d\boldsymbol{\mu}\_1 d\mathbf{y},\tag{28}$$

$$\mathcal{H}(t, \mathbf{x}) = \int\_{\mathbf{x}}^{\infty} \int\_{0}^{\mathbf{x}} \frac{\phi(\mathbf{u})}{\phi(\mathbf{y})} b(\mathbf{u}|\mathbf{y}) S(\mathbf{y}) n(t, \mathbf{y}) d\mathbf{u} d\mathbf{y}. \tag{29}$$

In a similar manner, we can represent the three-dimensional hypervolume conservation equations. Consider that *<sup>m</sup>*(*t*, **<sup>x</sup>**) := *<sup>ψ</sup>*(**x**)*f*(*t*, **<sup>x</sup>**) is the solution function, and the breakage function follows the relation (20). Then simply by replacing the function *φ* by *ψ* and simultaneously taking care of all the corresponding changes in the Equation (22), one can easily represent the three-dimensional hypervolume conservative model. Furthermore, one can extend the conservative formulations for problems with *n*-number of particle property components.

#### **3. Numerical Formulations**

#### *3.1. Two-Dimensional Model*

In this section, we present the discretized form of Equation (9). For this purpose, the truncated rectangular domain considered is *<sup>V</sup>* :=]0, *<sup>X</sup>*1]×]0, *<sup>X</sup>*2]. Let *<sup>I</sup>*<sup>1</sup> and *<sup>I</sup>*<sup>2</sup> be two positive integers, and *<sup>V</sup>* is further discretized in (*I*<sup>1</sup> <sup>×</sup> *<sup>I</sup>*2) number of rectangular subcells *<sup>V</sup>***<sup>i</sup>** :<sup>=</sup> \* *xi*1<sup>−</sup>1/2, *xi*<sup>1</sup>+1/2\* × \* *xi*2<sup>−</sup>1/2, *xi*<sup>2</sup>+1/2\* , where **<sup>i</sup>** :<sup>=</sup> (*i*1, *<sup>i</sup>*2), *<sup>I</sup>* :<sup>=</sup> (*I*1, *<sup>I</sup>*2) such that 1 <sup>≤</sup> **<sup>i</sup>** <sup>≤</sup> *<sup>I</sup>* along with *<sup>x</sup>*1/2 <sup>=</sup> *<sup>y</sup>*1/2 <sup>=</sup> 0, and *xI*<sup>1</sup>+1/2 <sup>=</sup> *<sup>X</sup>*1, *xI*<sup>2</sup>+1/2 <sup>=</sup> *<sup>X</sup>*2. Let <sup>Δ</sup>*i*<sup>1</sup> :<sup>=</sup> *xi*<sup>1</sup>+1/2 <sup>−</sup> *xi*1<sup>−</sup>1/2, <sup>Δ</sup>*i*<sup>2</sup> :<sup>=</sup> *xi*<sup>2</sup>+1/2 <sup>−</sup> *xi*2−1/2 and <sup>Δ</sup>**<sup>i</sup>** :<sup>=</sup> <sup>Δ</sup>*i*1Δ*i*<sup>2</sup> . Further, let **xi** :<sup>=</sup> *xi*<sup>1</sup> , *xi*<sup>2</sup> be the pivot or representative of the cell *V***i**, and the components of **xi** are defined by

$$\mathfrak{x}\_{i\_1} := \frac{\mathfrak{x}\_{i\_1 + 1/2} - \mathfrak{x}\_{i\_1 - 1/2}}{2}, \qquad \mathfrak{x}\_{i\_2} := \frac{\mathfrak{x}\_{i\_2 + 1/2} - \mathfrak{x}\_{i\_2 - 1/2}}{2}.$$

Under the above considerations, the flux flow at the right boundaries of the cell are given by F *xi*<sup>1</sup>+1/2, *xi*<sup>2</sup> , *<sup>t</sup>* , G *xi*<sup>1</sup> , *xi*<sup>2</sup>+1/2, *<sup>t</sup>* and H *xi*<sup>1</sup>+1/2, *xi*<sup>2</sup>+1/2, *<sup>t</sup>* , and similarly the flux flow at the other boundaries are defined.

Let *<sup>n</sup>***<sup>i</sup>** be the average value of the solution *<sup>n</sup>*(*t*, **<sup>x</sup>**) over the cell *<sup>V</sup>***i**, and is defined by

$$m\_{\mathbf{i}} = \frac{1}{\Delta\_{\mathbf{i}}} \int\_{V\_{\mathbf{i}}} n(t, \mathbf{x}) \, \mathbf{dx}. \tag{30}$$

.

Consider that *<sup>n</sup>*ˆ**i**(*t*) denotes the numerical approximation of *<sup>n</sup>***i**. For notational convenience, we drop the argument of the *<sup>t</sup>* from *<sup>n</sup>*ˆ**i**(*t*) in further discussions and simply denote it as *n*ˆ**i**.

Let us now evaluate the numerical approximation of the flux F *xi*<sup>1</sup>+1/2, *xi*<sup>2</sup> , *<sup>t</sup>* 

$$\begin{split} \mathcal{F}\left(\mathbf{x}\_{i\_{1}+1/2},\mathbf{x}\_{i\_{2}},t\right) &= \int\_{\mathbf{x}\_{i\_{1}+1/2}}^{\mathbf{x}\_{1}} \int\_{\mathbf{x}\_{i\_{2}}}^{\mathbf{x}\_{2}} \left[ \int\_{0}^{\mathbf{x}\_{i\_{1}+1/2}} \left(\boldsymbol{\epsilon} + \mathbf{x}\_{i\_{2}}\right) b\left(\boldsymbol{\epsilon},\mathbf{x}\_{i\_{2}}|\mathbf{u}\right) \mathbf{d}\boldsymbol{\epsilon} \right] \frac{S(\mathbf{u})}{\boldsymbol{\Phi}(\mathbf{u})} n(t,\mathbf{u}) d\mathbf{u} \\ &= \sum\_{k\_{1}=i\_{1}+1}^{I\_{1}} \int\_{\mathbf{x}\_{k\_{1}-1/2}}^{\mathbf{x}\_{k\_{1}+1/2}} \sum\_{k\_{2}=i\_{2}}^{I\_{2}} \int\_{\boldsymbol{\beta}\left(i\_{2}k\_{2}\right)}^{\mathbf{x}\_{k\_{2}+1/2}} \left[ \sum\_{l\_{1}=1}^{i\_{1}} \int\_{\mathbf{x}\_{l\_{1}-1/2}}^{\mathbf{x}\_{l\_{1}+1/2}} \left(\boldsymbol{\epsilon} + \mathbf{x}\_{i\_{2}}\right) b\left(\boldsymbol{\epsilon},\mathbf{x}\_{i\_{2}}|\mathbf{u}\right) d\mathbf{d}\boldsymbol{\epsilon} \right] \\ & \qquad \times \frac{S(\mathbf{u})}{\boldsymbol{\Phi}(\mathbf{u})} n(t,\mathbf{u}) d\mathbf{u}, \end{split}$$

Here,

$$\beta(i\_2, k\_2) := \begin{cases} \ x\_{k\_{2'}} & \text{when} \quad i\_2 = k\_{2'}, \\\ x\_{k\_2 - 1/2} & \text{otherwise}. \end{cases}$$

Applying quadrature formulae to the integrals, the numerical flux is given by

$$\mathcal{F}\_{i\_1+1/2,i\_2} := \sum\_{\substack{l\_1=i\_1+1 \ k\_2=i\_2}}^{I\_1} \sum\_{k\_2=i\_2}^{I\_2} \hat{m}\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^{\mathcal{S}} \sum\_{l\_1=1}^{i\_1} \mathcal{B}\_{l\_1,i\_2|\mathbf{k}} \Delta\_{l\_1} \Delta\_{\mathbf{k}\_{\mathbf{k}}} \tag{31}$$

where **<sup>k</sup>** :<sup>=</sup> (*k*1, *<sup>k</sup>*2), **<sup>l</sup>** :<sup>=</sup> (*l*1, *<sup>l</sup>*2),

$$\mathcal{B}\_{l\_1,j\_2|\mathbf{k}} := (\mathbf{x}\_{l\_1} + y\_{i\_2})b\left(\mathbf{x}\_{l\_1}, y\_{i\_2}|\mathbf{x}\_{\mathbf{k}}\right), \quad \text{and} \quad \mathcal{A}\_{\mathbf{k}}^{\mathcal{B}} := \int\_{\mathbf{x}\_{k\_1-1/2}}^{\mathbf{x}\_{k\_1+1/2}} \int\_{\beta(i\_2k\_2)}^{y\_{k\_2+1/2}} \frac{\mathcal{S}(\mathbf{u})}{\phi(\mathbf{u})} d\mathbf{u}.$$

In a similar manner, under the following notations

$$\mathcal{B}\_{i\_1,l\_2|\mathbf{k}} := \left(\mathbf{x}\_{l\_1} + \mathbf{x}\_{l\_2}\right) b\left(\mathbf{x}\_{l\_1}, \mathbf{x}\_{l\_2}|\mathbf{x}\_{\mathbf{k}}\right), \quad \mathcal{B}\_{l|\mathbf{k}} := \left(\mathbf{x}\_{l\_1} + \mathbf{x}\_{l\_2}\right) b\left(\mathbf{x}\_{l\_1}, y\_{l\_2}|\mathbf{x}\_{\mathbf{k}}\right).$$

$$\mathcal{A}\_{\mathbf{k}}^{a} := \int\_{a(i\_1k\_1)}^{x\_{k\_1+1/2}} \int\_{x\_{k\_2-1/2}}^{x\_{k\_2+1/2}} \frac{\mathcal{S}(\mathbf{u})}{\varPhi(\mathbf{u})} d\mathbf{u}, \quad \text{and} \quad \mathcal{A}\_{\mathbf{k}} := \int\_{x\_{k\_1-1/2}}^{x\_{k\_1+1/2}} \int\_{\varprojlim\_{t \in \mathbf{u}} -1/2}^{\mathfrak{Y}k\_2 + 1/2} \frac{\mathcal{S}(\mathbf{u})}{\varPhi(\mathbf{u})} d\mathbf{u}, \quad \text{and} \quad \mathfrak{X} := \int\_{\varprojlim\_{t \in \mathbf{u}} -1/2}^{\mathfrak{Y}k\_2 + 1/2} \frac{\mathcal{S}(\mathbf{u})}{\varprojlim\_{t \in \mathbf{u}} -1/2} d\mathbf{u}.$$

with

$$\mathfrak{a}(i\_1, k\_1) := \begin{cases} \begin{array}{c} \mathfrak{x}\_{k\_1 \prime} \\ \updownarrow\_{k\_1 - 1/2\prime} \end{array} & \text{when } i\_1 = k\_1 \text{.} \\\end{cases}$$

the other numerical fluxes at the cell interfaces are written as

$$\mathcal{G}\_{i\_1, i\_2+1/2} := \sum\_{k\_1=i\_1}^{I\_1} \sum\_{k\_2=i\_2+1}^{I\_2} \hat{n}\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^a \sum\_{l\_2=1}^{i\_2} \mathcal{B}\_{i\_1 J\_2 \vert \mathbf{k}} \Delta\_{l\_2} \Delta\_{\mathbf{k} \prime} \tag{32}$$

and

$$\mathcal{H}\_{\mathbf{i}+1/2} := \sum\_{\mathbf{k}=\mathbf{i}+1}^{l} \boldsymbol{\hat{n}}\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{l}=1}^{\mathbf{i}} \mathcal{B}\_{\mathbf{l}|\mathbf{k}} \Delta\_{\mathbf{l}} \Delta\_{\mathbf{k}}.\tag{33}$$

Therefore, the semi-discrete finite volume representation of Equation (9) is written as

$$\begin{split} \frac{d\boldsymbol{\hbar}\_{1}}{dt}\Delta\_{\mathrm{i}} &= \left[\mathcal{F}\_{i\_{1}+1/2,i\_{2}} - \mathcal{F}\_{i\_{1}-1/2,i\_{2}}\right] \Delta\_{i\_{2}} + \left[\mathcal{G}\_{i\_{1},i\_{2}+1/2} - \mathcal{G}\_{i\_{1},i\_{2}-1/2}\right] \Delta\_{i\_{1}} \\ &- \left[\mathcal{H}\_{i\_{1}+1/2,i\_{2}+1/2} - \mathcal{H}\_{i\_{1}+1/2,i\_{2}-1/2} - \mathcal{H}\_{i\_{1}-1/2,i\_{2}+1/2} + \mathcal{H}\_{i\_{1}-1/2,i\_{2}-1/2}\right]. \end{split} \tag{34}$$

The above scheme (34) obeys the discrete volume conservation law (the detailed calculations are given in Appendix A). However, in the subsequent section, we shall numerically validate that scheme (34) fails to predict the evolution of total number of fragments with good accuracy. In this context, the flux functions are redefined by introducing a weight function which enables the model to obey volume conservation laws, as well as predict the zeroth moment with high accuracy. The newly proposed semi-discrete formulation is written as follows:

$$\begin{split} \frac{d\hat{\mathbf{n}}\_{1}}{dt}\Delta\_{1} &= \left[\hat{\mathcal{F}}\_{i\_{1}+1/2,i\_{2}} - \hat{\mathcal{F}}\_{i\_{1}-1/2,i\_{2}}\right]\Delta\_{i\_{2}} + \left[\hat{\mathcal{G}}\_{i\_{1}i\_{2}+1/2} - \hat{\mathcal{G}}\_{i\_{1}i\_{2}-1/2}\right]\Delta\_{i\_{1}} \\ &- \left[\hat{\mathcal{H}}\_{i\_{1}+1/2,i\_{2}+1/2} - \hat{\mathcal{H}}\_{i\_{1}+1/2,i\_{2}-1/2} - \hat{\mathcal{H}}\_{i\_{1}-1/2,i\_{2}+1/2} + \hat{\mathcal{H}}\_{i\_{1}-1/2,i\_{2}-1/2}\right], \end{split} \tag{35}$$

where the modified fluxes at the cell boundaries are defined as follows:

$$\hat{\mathcal{F}}\_{i\_1+1/2,i\_2} := \sum\_{k\_1=i\_1+1}^{I\_1} \sum\_{k\_2=i\_2}^{I\_2} \hat{m}\_{\mathbf{k}} \delta\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^{\S} \sum\_{l\_1=1}^{i\_1} \mathcal{B}\_{l\_1,i\_2|\mathbf{k}} \Delta\_{l\_1} \Delta\_{\mathbf{k}\prime} \tag{36}$$

$$\mathcal{G}\_{i\_1, i\_2+1/2} := \sum\_{k\_1=i\_1}^{I\_1} \sum\_{k\_2=i\_2+1}^{I\_2} \hat{m}\_{\mathbf{k}} \delta\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^a \sum\_{l\_2=1}^{i\_2} \mathcal{B}\_{i\_1, l\_2 | \mathbf{k}} \Delta\_{l\_2} \Delta\_{\mathbf{k}'} \tag{37}$$

$$\mathcal{H}\_{\mathbf{i}+1/2} := \sum\_{\mathbf{k}=\mathbf{i}+1}^{l} \boldsymbol{\hbar}\_{\mathbf{k}} \boldsymbol{\delta}\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{l}=1}^{\mathbf{i}} \mathcal{B}\_{\mathbf{l}|\mathbf{k}} \Delta\_{\mathbf{l}} \Delta\_{\mathbf{k}\prime} \tag{38}$$

and *δ* is the weight factor, defined by

$$\delta\_{\mathbf{k}} := \frac{S\_{\mathbf{k}}[\boldsymbol{\nu}(\mathbf{x\_{k}}) - 1]}{\mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{l}=1}^{\mathbf{k}} (\boldsymbol{\phi}(\mathbf{x\_{k}}) - \boldsymbol{\phi}(\mathbf{x\_{l}})) \mathcal{B}\_{\mathbf{i}|\mathbf{k}}},\tag{39}$$

along with *<sup>δ</sup>*1,1 <sup>=</sup> 1. In the above definition of the weight, the terms *<sup>S</sup>***<sup>k</sup>** and *<sup>ν</sup>*(**xk**) denote the discrete selection function and the number of fragments, respectively.

In the following section, we will numerically validate that the two-dimensional scheme (35) is consistent with the zeroth moment and it also obeys the volume conservation law. The proof of this claim follows similar to that of the model (34) and can easily be followed from the outline given in Appendix A.

In a similar manner, a new scheme preserving the cluster hypervolume and estimating the continuous model (16) along with the zeroth moment can be defined as follows:

$$\begin{split} \frac{d\mathbf{\hat{m}}\_{1}}{dt}\Delta\_{\mathbf{i}} &= \left[\hat{\mathcal{F}}\_{i\_{1}+1/2,j\_{2}} - \hat{\mathcal{F}}\_{i\_{1}-1/2,j\_{2}}\right] \Delta\_{i\_{2}} + \left[\hat{\mathcal{G}}\_{i\_{1},j\_{2}+1/2} - \hat{\mathcal{G}}\_{i\_{1},j\_{2}-1/2}\right] \Delta\_{i\_{1}} \\ &- \left[\hat{\mathcal{H}}\_{i\_{1}+1/2,j\_{2}+1/2} - \hat{\mathcal{H}}\_{i\_{1}+1/2,j\_{2}-1/2} - \hat{\mathcal{H}}\_{i\_{1}-1/2,j\_{2}+1/2} + \hat{\mathcal{H}}\_{i\_{1}-1/2,j\_{2}-1/2}\right], \end{split} \tag{40}$$

where the discrete flux functions at the cell boundaries are defined by

$$\hat{\mathcal{F}}\_{i\_1+1/2,i\_2} := \sum\_{\substack{k\_1=i\_1+1 \ k\_2=i\_2}}^{I\_1} \sum\_{k\_2=i\_2}^{I\_2} \eta \mathbf{\hat{h}}\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^{\S} \sum\_{l\_1=1}^{i\_1} \mathcal{B}\_{l\_1,i\_2|\mathbf{k}} \Delta\_{l\_1} \Delta\_{\mathbf{k}\_{\mathbf{k}}} \tag{41}$$

$$\hat{\mathcal{G}}\_{i\_1, i\_2+1/2} := \sum\_{k\_1=i\_1}^{I\_1} \sum\_{k\_2=i\_2+1}^{I\_2} \hat{m}\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^{\underline{a}} \sum\_{l\_2=1}^{i\_2} \mathcal{B}\_{i\_1, l\_2 \mid \mathbf{k}} \Delta\_{l\_2} \Delta\_{\mathbf{k} \prime} \tag{42}$$

$$\hat{\mathcal{H}}\_{\mathbf{i}+1/2} := \sum\_{\mathbf{k}=\mathbf{i}+1}^{I} \hat{m}\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{l}=1}^{\mathbf{i}} \mathcal{B}\_{\mathbf{l}|\mathbf{k}} \Delta\_{\mathbf{l}} \Delta\_{\mathbf{k}\mathbf{l}} \tag{43}$$

with the discrete breakage function as

$$\mathcal{B}\_{\mathbf{i}|\mathbf{k}} := \psi(\mathbf{x\_i}) b(\mathbf{x\_i}|\mathbf{x\_k})\_{\prime}$$

and *ω* is the weight factor, defined by

$$
\omega\_{\mathbf{k}} := \frac{S\_{\mathbf{k}}[\nu(\mathbf{x\_k}) - 1]}{\mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{i}=1}^{\mathbf{k}} (\psi(\mathbf{x\_k}) - \psi(\mathbf{x\_i})) \mathcal{B}\_{\mathbf{i}|\mathbf{k}}},\tag{44}
$$

**Remark 1.** *In the proposed two-dimensional model* (10)*, there are three numerical fluxes operating at the cell boundaries. An interesting feature of the new two-dimensional model is that a single weight function is sufficient for redefining the modified scheme to become consistent with the zerothand the first-order moments.*

**Remark 2.** *It is to be noted that the finite volume scheme* (34) *is obtained by direct application of the midpoint quadrature rules to the continuous Equation* (9)*. Thus, it represents the numerical model of [23] with two degrees of freedom.*

### *3.2. Three-Dimensional Model*

*r*=1

In this part, we present the three-dimensional finite volume scheme approximating the multi-fragmentation model. Here, the scheme is expressed using vector notation, which will give an outline for further extension of the proposed scheme in higher dimensions.

Similar to the two-dimensional model, the computational domain considered is *<sup>V</sup>* := 3 ∏]0, *Xr*] which is further divided into a finite number of sub-cells

$$V\_{\mathbf{i}} := \prod\_{r=1}^{3} [\mathfrak{x}\_{i\_r - 1/2, r} \mathfrak{x}\_{i\_r + 1/2}]$$

with *ir* <sup>=</sup> 1, 2, ... , *Ir*. Let *<sup>n</sup>*ˆ**<sup>i</sup>** be the numerical approximation of *<sup>n</sup>*(*t*, **<sup>x</sup>**) over the cell *<sup>V</sup>***i**. Further, let Δ**<sup>i</sup>** denote the volume of the cell *V***i**, and the the cell representative is defined by **xi** :<sup>=</sup> *xi*<sup>1</sup> , *xi*<sup>2</sup> , *xi*<sup>3</sup> . Consider that the breakage function obeys the conservative Formula (14), then the three-dimensional extension of the newly proposed volumeconservative formulation (33) is written as

d*n*ˆ**<sup>i</sup>** d*t* <sup>Δ</sup>**<sup>i</sup>** <sup>=</sup> <sup>F</sup><sup>ˆ</sup> (1) *<sup>i</sup>*1+1/2,*i*2,*i*<sup>3</sup> <sup>−</sup> <sup>F</sup><sup>ˆ</sup> (1) *i*1−1/2,*i*2,*i*<sup>3</sup> <sup>Δ</sup>*i*2,*i*<sup>3</sup> <sup>+</sup> <sup>F</sup><sup>ˆ</sup> (2) *<sup>i</sup>*1,*i*2+1/2,*i*<sup>3</sup> <sup>−</sup> <sup>F</sup><sup>ˆ</sup> (2) *i*1,*i*2−1/2,*i*<sup>3</sup> Δ*i*1,*i*<sup>3</sup> + <sup>F</sup><sup>ˆ</sup> (3) *<sup>i</sup>*1,*i*2,*i*3+1/2 <sup>−</sup> <sup>F</sup><sup>ˆ</sup> <sup>3</sup> *<sup>i</sup>*1,*i*2,*i*3−1/2 Δ*i*1,*i*<sup>2</sup> − Gˆ(1) *<sup>i</sup>*1,*i*2+1/2,*i*3+1/2 <sup>−</sup> <sup>G</sup>ˆ(1) *<sup>i</sup>*1,*i*2−1/2,*i*3+1/2 <sup>−</sup> <sup>G</sup>ˆ(1) *<sup>i</sup>*1,*i*2+1/2,*i*3−1/2 <sup>+</sup> <sup>G</sup>ˆ(1) *<sup>i</sup>*1,*i*2−1/2,*i*3−1/2 Δ*i*1 − Gˆ(2) *<sup>i</sup>*1+1/2,*i*2,*i*3+1/2 <sup>−</sup> <sup>G</sup>ˆ(2) *<sup>i</sup>*1−1/2,*i*2,*i*3+1/2 <sup>−</sup> <sup>G</sup>ˆ(2) *<sup>i</sup>*1+1/2,*i*2,*i*3−1/2 <sup>+</sup> <sup>G</sup>ˆ(2) *<sup>i</sup>*1−1/2,*i*2,*i*3−1/2 Δ*i*2 − Gˆ(3) *<sup>i</sup>*1+1/2,*i*2+1/2,*i*<sup>3</sup> <sup>−</sup> <sup>G</sup>ˆ(3) *<sup>i</sup>*1−1/2,*i*2+1/2,*i*<sup>3</sup> <sup>−</sup> <sup>G</sup>ˆ(3) *<sup>i</sup>*1+1/2,*i*2−1/2,*i*<sup>3</sup> <sup>+</sup> <sup>G</sup>ˆ(3) *i*1−1/2,*i*2−1/2,*i*<sup>3</sup> Δ*i*3 + ) <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1+1/2,*i*2+1/2,*i*3+1/2 <sup>−</sup> <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1+1/2,*i*2+1/2,*i*3−1/2 <sup>−</sup> <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1+1/2,*i*2−1/2,*i*3+1/2 <sup>+</sup> <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1+1/2,*i*2−1/2,*i*3−1/2 <sup>−</sup> <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1−1/2,*i*2+1/2,*i*3+1/2 <sup>+</sup> <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1−1/2,*i*2−1/2,*i*3+1/2 <sup>+</sup>H<sup>ˆ</sup> *<sup>i</sup>*1−1/2,*i*2+1/2,*i*3−1/2 <sup>−</sup> <sup>H</sup><sup>ˆ</sup> *<sup>i</sup>*1−1/2,*i*2−1/2,*i*3−1/2\* . (45)

Considering **<sup>k</sup>** :<sup>=</sup> (*k*1, *<sup>k</sup>*2, *<sup>k</sup>*3) and **<sup>l</sup>** :<sup>=</sup> (*l*1, *<sup>l</sup>*2, *<sup>l</sup>*3), the redefined flux functions are written as follows:

$$\mathcal{F}^{(1)}\_{i\_1+1/2,i\_2;j\_3} := \sum\_{k\_1=i\_1+1}^{I\_1} \sum\_{\substack{(k\_2,k\_3)=(i\_2,i\_3)}}^{(I\_2,I\_3)} \hat{n}\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}^{\S,\gamma}\_{\mathbf{k}} \Delta\_{\mathbf{k}} \sum\_{l\_1=1}^{i\_1} \mathcal{B}\_{l\_1,i\_2,i\_3|\mathbf{k}} \Delta\_{l\_1} \tag{46}$$

$$\mathcal{F}\_{i\_1,i\_2+1/2,i\_3}^{(2)} := \sum\_{\substack{k\_2=i\_2+1 \ (k\_1,k\_3)=(i\_1,i\_3)}}^{l\_2} \sum\_{\substack{(k\_1,i\_3)=(i\_1,i\_3)}}^{I\_{1,I\_3}} \beta\_{\mathbf{k}} \omega\_{\mathbf{k}'} \mathcal{A}\_{\mathbf{k}}^{a,\gamma} \Delta\_{\mathbf{k}} \sum\_{l\_2=1}^{i\_2} \mathcal{B}\_{i\_1,l\_2,i\_3|\mathbf{k}} \Delta\_{l\_2 \prime} \tag{47}$$

$$\mathcal{F}\_{i\_1,i\_2,i\_3+1/2}^{(3)} := \sum\_{(k\_1,k\_2)=(i\_1,i\_2)}^{(I\_1,I\_2)} \sum\_{k\_3=i\_3+1}^{I\_3} \hbar\_{\mathbf{k}} \omega\_{\mathbf{k}'} \mathcal{A}\_{\mathbf{k}}^{a,\emptyset} \Delta\_{\mathbf{k}} \sum\_{l\_3=1}^{i\_3} \mathcal{B}\_{i\_1,i\_2,l\_3|\mathbf{k}} \Delta\_{l\_3 \prime} \tag{48}$$

$$\mathcal{G}^{(1)}\_{i\_1j\_2+1/2,j\_3+1/2} := \sum\_{k\_1=i\_1}^{I\_1} \sum\_{\substack{(k\_2,k\_3)=(i\_2+1,i\_3+1)}}^{(I\_2,l\_3)} \hat{n}\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}^{\mathbf{a}}\_{\mathbf{k}} \Delta\_{\mathbf{k}} \sum\_{\substack{(l\_2,l\_3)=1}}^{i\_2,i\_3} \mathcal{B}\_{i\_1l\_2l\_3|\mathbf{k}} \Delta\_{l\_2l\_3r} \tag{49}$$

$$\mathcal{G}^{(2)}\_{i\_1+1/2,i\_2,i\_3+1/2} := \sum\_{(k\_1,k\_3=(i\_1+1,i\_3+1))}^{(l\_1,l\_3)} \sum\_{k\_2=i\_2}^{l\_2} \, \hbar\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}^{\emptyset}\_{\mathbf{k}} \Delta\_{\mathbf{k}} \sum\_{(l\_1,l\_2)=1}^{i\_1,i\_3} \mathcal{B}\_{l\_1,l\_2,i\_3|\mathbf{k}} \Delta\_{l\_1,l\_2,i} \tag{50}$$

$$\hat{\mathcal{G}}\_{i\_1+1/2,i\_2+1/2,j\_3}^{(3)} := \sum\_{\substack{k\_1=i\_1+1 \ (k\_2,k\_3)=(i\_2+1,j\_3+1)}}^{l\_1} \hat{\eta}\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}}^{\gamma} \Delta\_{\mathbf{k}} \sum\_{(l\_1,l\_2)=1}^{i\_1j\_2} \mathcal{B}\_{l\_1,l\_2,j\_3|\mathbf{k}} \Delta\_{l\_2,l\_3},\tag{51}$$

$$
\hat{\mathcal{H}}\_{\mathbf{i}+1/2} := \sum\_{\mathbf{k}=\mathbf{i}}^{\mathrm{I}} \hbar\_{\mathbf{k}} \omega\_{\mathbf{k}} \mathcal{A}\_{\mathbf{k}} \Delta\_{\mathbf{k}} \sum\_{\mathbf{l}=\mathbf{l}}^{\mathrm{I}} \mathcal{B}\_{\mathbf{l}|\mathbf{k}} \Delta\_{\mathbf{l}}.\tag{52}
$$

Here, *δ***<sup>k</sup>** is the weight factor, defined by

$$\delta\_{\mathbf{k}} := \frac{S\_{\mathbf{k}}[\boldsymbol{\nu}(\mathbf{x\_{k}}) - 1]}{\mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{l}=1}^{\mathbf{k}} [\boldsymbol{\phi}(\mathbf{x\_{k}}) - \boldsymbol{\phi}(\mathbf{x\_{l}})] \mathcal{B}\_{\mathbf{l}|\mathbf{k}}},\tag{53}$$

and the factors <sup>A</sup>*α*,*<sup>β</sup>* **<sup>k</sup>** , <sup>A</sup>*β*,*<sup>γ</sup>* **<sup>k</sup>** and <sup>A</sup>*α*,*<sup>γ</sup>* **<sup>k</sup>** are defined as follows:

$$\mathcal{A}\_{\mathbf{k}}^{\alpha,\emptyset} := \int\_{\alpha(i\_1,l\_1)}^{\chi\_{l\_1+1/2}} \int\_{\beta(i\_2,l\_2)}^{\chi\_{l\_2+1/2}} \int\_{x\_{l\_3-1/2}}^{\chi\_{l\_3+1/2}} \frac{S(\mathbf{u})}{\phi(\mathbf{u})} \mathbf{d} \mathbf{u} \,\mathrm{d}\mathbf{u}$$

$$\begin{aligned} \mathcal{A}^{\S,\gamma}\_{\mathbf{k}} &:= \int\_{\chi\_{l\_1-1/2}}^{\chi\_{l\_1+1/2}} \int\_{\beta(i\_2,l\_2)}^{\chi\_{l\_2+1/2}} \int\_{\gamma(i\_3,l\_3)}^{\chi\_{l\_3+1/2}} \frac{S(\mathbf{u})}{\phi(\mathbf{u})} d\mathbf{u} \,\mathrm{d}\mathbf{u} \\\\ \mathcal{A}^{\mathbf{a},\gamma}\_{\mathbf{k}} &:= \int\_{\kappa(i\_1,l\_1)}^{\chi\_{l\_1+1/2}} \int\_{\chi\_{l\_2-1/2}}^{\chi\_{l\_2+1/2}} \int\_{\gamma(i\_3,l\_3)}^{\chi\_{l\_3+1/2}} \frac{S(\mathbf{u})}{\phi(\mathbf{u})} d\mathbf{u} \,\mathrm{d}\mathbf{u} \end{aligned}$$

along with <sup>A</sup>*<sup>α</sup>* **<sup>k</sup>**, <sup>A</sup>*<sup>β</sup>* **<sup>k</sup>**, <sup>A</sup>*<sup>γ</sup>* **<sup>k</sup>**, <sup>A</sup>**<sup>k</sup>** and *<sup>α</sup>*(*i*1, *<sup>l</sup>*1), *<sup>β</sup>*(*i*2, *<sup>l</sup>*2), *<sup>γ</sup>*(*i*3, *<sup>l</sup>*3) being defined in a similar manner as done before.

Following the same trail, one can easily define the three-dimensional hypervolume preserving numerical model. In this case, the numerical solution will be given by *m*ˆ **<sup>i</sup>**, which is the approximation of *<sup>m</sup>*(*t*, **<sup>x</sup>**) over the cell *<sup>V</sup>***<sup>i</sup>** and the breakage function will obey the hypervolume conservation law (20). Thus, the weight function will be defined as

$$
\omega\_{\mathbf{k}} := \frac{S\_{\mathbf{k}}[\nu(\mathbf{x\_k}) - 1]}{\mathcal{A}\_{\mathbf{k}} \sum\_{\mathbf{i}=1}^{\mathbf{k}} [\psi(\mathbf{x\_k}) - \psi(\mathbf{x\_i})] \mathcal{B}\_{\mathbf{i}|\mathbf{k}}}.\tag{54}
$$

#### **4. Results**

In this section, we validate the efficiency of the newly proposed finite volume models with the standard finite volume scheme over several test problems. Since the new models are defined with the help of a weight factor, we call it weighted finite volume scheme (WFVS). On the other hand, Remark 2 indicates that the standard forms of the schemes which are directly derived from the continuous equations were initially proposed by [23] for fragmentation models with one degree of freedom. Therefore, for future reference, we call the standard models the existing finite volume schemes (EFVS). However, we need to mention that the two-dimensional extension of EFVS is not available in the literature to date, and this article not only proposes an improved model, but also extends the existing finite volume schemes for multidimensional fragmentation events.

For one-dimensional fragmentation problems, Ref. [30,31] has obtained the exact solutions for a certain class of fragmentation kinetics. However, exact solutions in closed forms are very rare in the multidimensional setup. In order to validate the accuracy of the proposed schemes, we choose four test problems with two degrees of freedom, and two test problems with three degrees of freedom. For all the test problems, exact solutions in closed form are not always available in the literature. However, the zeroth and the first-order moments can be computed exactly, which is sufficient to validate the accuracy of the new schemes. Therefore, in the following section we discuss the efficiency of the WFV scheme to predict the different physically important moment functions over the EFV scheme. Our study builds up on both qualitative and quantitative assessments. The qualitative accuracy is represented through graphical representation of the different entities whereas, the qualitative analysis is performed by computing relative errors of the moment functions over different grid points.

The computational domain *<sup>V</sup>* <sup>=</sup> [0, 1] <sup>×</sup> [0, 1] is considered for all the two-dimensional test problems. The domain *V* is discretized into 15 × 15 nonuniform subintervals bearing the geometric relation *xi*+1/2 <sup>=</sup> *rxi*−1/2 where *<sup>r</sup>* :<sup>=</sup> 3.9811 is the geometric ratio. Additionally, all the test problems are supported by a mono-dispersed initial condition

$$f(\mathbf{x}\_1, \mathbf{x}\_2, 0) = \delta(\mathbf{x}\_1 - 1)\delta(\mathbf{x}\_2 - 1).$$

Similar extensions of the above data are considered for the three dimensional models. The semi-discrete schemes (7) and (8) are solved using MATLAB-2019B software in a standard PC with i5-7500 CPU processor @ 3.41 GHz and 8 GB RAM.

#### *4.1. Examples in Two-Dimensions*

#### 4.1.1. Volume Conservation Problems

In the first instance, we consider test problems with constant particle selection rate, that is *<sup>S</sup>*(*x*, *<sup>y</sup>*) = 1 and two different daughter distribution functions

$$b\_1(\mathbf{x}, y|\mathbf{x}', y') = \frac{2}{\mathbf{x}'y'} \quad \text{and} \quad b\_2\left(\mathbf{x}, y|\mathbf{x}', y'\right) = 2\delta\left(\mathbf{x} - \frac{\mathbf{x}'}{2}\right)\delta\left(y - \frac{y'}{2}\right).$$

The first breakage function *b*<sup>1</sup> is a size-independent function of its arguments and physically represents random scission of particles. On the other hand, the second breakage function *b*<sup>2</sup> represents size-dependent distribution of daughter fragments, choosing the daughter-particles exactly half the size of the parent particle. The exact solution in closed form for the above set of fragmentation kernels are not available in the literature. However, we can calculate the zeroth <sup>M</sup>0,0(*t*), first <sup>M</sup>1,0(*t*) +M0,1(*t*) and the cross moments <sup>M</sup>1,1(*t*) exactly for both the problems (calculations of the exact moments are given in Appendix B). In this regard, the exact moments are given in the Table 1.

**Table 1.** Exact moment functions for *<sup>S</sup>* = 1 and breakage functions *<sup>b</sup>*1, *<sup>b</sup>*2.


Figures 1 and 2 represent the numerical moments obtained from EFVS and WFVS against their exact values with breakage functions *b*<sup>1</sup> and *b*2, respectively. More precisely, Figures 1a and 2a present the comparison of zeroth and first-order moment functions, and Figures 1b and 2b presents the first-order cross moment function <sup>M</sup>1,1(*t*). In order to obtain a clear visibility of different markers, we plot Figures 1a and 2a on a semilogarithmic scale with respect to the *y*-axis. In both the cases, it is observed that WFVS estimates the zeroth-order, first-order and the cross moments with high accuracy, whereas the EFVS conserves the total volume but fails to produce a good estimate of the other moments.

(**a**) Zeroth- and first-order moments.

**Figure 1.** *Cont*.

(**b**) First-order cross moments.

**Figure 1.** Comparison of different moments with selection function *<sup>S</sup>* <sup>=</sup> 1, and breakage function *<sup>b</sup>*1.

(**b**) First-order cross moments

**Figure 2.** Comparison of different moments with selection function *<sup>S</sup>* = 1, and breakage function *<sup>b</sup>*2.

In the Table 2, the relative error of the moment functions for *<sup>S</sup>* <sup>=</sup> 1 and *<sup>b</sup>*<sup>1</sup> are calculated at *<sup>t</sup>* = 10 over three different grid sizes. The geometric ratios to generate 10 <sup>×</sup> 10, 15 <sup>×</sup> <sup>15</sup> and 20 <sup>×</sup> 20 grids are 7.9433, 3.9811 and 2.8184, respectively. The discrete *<sup>L</sup>*1-error norm

$$\text{error} := \sum\_{i=1}^{I} \left| \frac{\mathcal{M}^{exact} - \mathcal{M}^{num}}{\mathcal{M}^{exact}} \right|.$$

is used to calculate the errors. Similarly, the relative error acquired while computing the moments for *<sup>S</sup>* <sup>=</sup> 1 and *<sup>b</sup>*<sup>2</sup> over different grid points are represented in Table 3.

**Table 2.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = 1 and breakage functions *<sup>b</sup>*<sup>1</sup> at *<sup>t</sup>* <sup>=</sup> 10.


**Table 3.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = 1 and breakage functions *<sup>b</sup>*<sup>2</sup> at *<sup>t</sup>* <sup>=</sup> 10.


In the second instance, we choose two problems by setting the size-dependent selection function *<sup>S</sup>*(*x*, *<sup>y</sup>*) = *<sup>x</sup>* + *<sup>y</sup>*, along with the previously chosen particle daughter distribution functions *b*<sup>1</sup> and *b*2. In this case, also the exact solutions are not available in the literature, however only the zeroth- and first-order moment functions can be calculated exactly. In both the cases, the moment functions are given as

$$
\mathcal{M}\_{1,0}(t) + \mathcal{M}\_{0,1}(t) = 1,\quad \text{and}\quad \mathcal{M}\_{0,0}(t) = 1 + 2t.
$$

The following Figure 3a,b represent the efficiency of the WFV scheme over the EFV scheme to estimate the zeroth- and the first-order moments. It is observed that both the schemes obey the volume conservation laws with high accuracy, but the WFV scheme is highly accurate to predict the evolution of total number of particles. Tables 4 and 5 represent the relative errors over different grid points at time *<sup>t</sup>* = 10.

**Figure 3.** Comparison of zeroth- and first-order moments with selection function *<sup>S</sup>* = *<sup>x</sup>* + *<sup>y</sup>*, and breakage functions *b*1, *b*2.

**Table 4.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* <sup>=</sup> *<sup>x</sup>* <sup>+</sup> *<sup>y</sup>* and breakage function *<sup>b</sup>*<sup>1</sup> at *<sup>t</sup>* <sup>=</sup> 10.


**Table 5.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = *<sup>x</sup>* + *<sup>y</sup>* and breakage function *<sup>b</sup>*<sup>2</sup> at *<sup>t</sup>* <sup>=</sup> 10.


### 4.1.2. Hypervolume Conservation Problems

In this part, the considered breakage functions *b* should satisfy the hypervolume conservation rule (8). In this regard, we choose the following breakage functions

$$b\_3(\mathbf{x}, y | \mathbf{x}', y') = \frac{4}{\mathbf{x}' y'} \quad \text{and} \quad b\_4(\mathbf{x}, y | \mathbf{x}', y') = \frac{\mathbf{x}' \delta(\mathbf{x} - \mathbf{x}') + y' \delta(y - y')}{\mathbf{x}' y'}.\tag{55}$$

The breakage function *b*<sup>3</sup> is independent of the daughter-particle size, whereas *b*<sup>4</sup> represents the particle breakage along the longer side of the rectangular structure. Similar to the examples of volume conservation models, we choose two types of selection functions: (i) size-independent kernels *<sup>S</sup>*(*x*, *<sup>y</sup>*) = 1, and (ii) size-dependent kernels *<sup>S</sup>*(*x*, *<sup>y</sup>*) = *xy*.

Like before, the exact solutions are not available in the literature, however we can calculate three moment functions exactly, and they are given in Table 6.

**Table 6.** Exact moment functions for *<sup>S</sup>* = 1 and breakage functions *<sup>b</sup>*3, *<sup>b</sup>*4.


In Figure 4a, we plot the zeroth- and the cross-moment functions and observe that the new WFV scheme predicts the corresponding moments with high accuracy. On the other hand, we plot the first-order moment in Figure 4b. In this case also, the weighted scheme exhibits high accuracy to estimate the moment compared to the standard model. In Figure 4b, we take the the axes in loglog scale for a distinct visibility of the plots. In this problem, Figure 4c represents the comparison of hypervolume distribution functions with the numerical values obtained from the two schemes. We follow the flat pictorial representation to plot the hypervolume distribution as presented in [32]. For the other problems, only the exact moment functions can be calculated for comparison with the numerical models.

The relative errors over different grid points are presented at Table 7.

In the second instance, we consider the size-dependent selection function *<sup>S</sup>*(*x*, *<sup>y</sup>*) = <sup>1</sup> and *b*<sup>4</sup> as the daughter distribution function. The exact solution is not available for this problem, but we can evaluate the zeroth and cross moments exactly. From the Figure 5 and Table 8, we can see that the newly proposed WFV scheme predicts the moments with high accuracy.

**Table 7.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = 1 and breakage functions *<sup>b</sup>*<sup>3</sup> at *<sup>t</sup>* <sup>=</sup> 10.


(**a**) Zeroth- and first-order cross moment.

(**c**) Particle size distribution.

**Figure 4.** Comparison of different moments with selection function *<sup>S</sup>* <sup>=</sup> 1, and breakage function *<sup>b</sup>*3.

(**a**) Zeroth- and first-order cross moment

**Figure 5.** Comparison of different moments with selection function *<sup>S</sup>* <sup>=</sup> 1, and breakage function *<sup>b</sup>*4.

**Table 8.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = 1 and breakage functions *<sup>b</sup>*<sup>4</sup> at *<sup>t</sup>* <sup>=</sup> 10.


Next, we consider two problems with size-dependent selection function *<sup>S</sup>*(*x*, *<sup>y</sup>*) = *xy* and the breakage functions *b*<sup>3</sup> and *b*4. Only the zeroth- and first-order moment functions can be calculated in closed form, and are given in the Table 9.

**Table 9.** Exact moment functions for *<sup>S</sup>* = *xy* and breakage functions *<sup>b</sup>*3, *<sup>b</sup>*4.


Numerical evaluation of the moments using the WFV and EFV schemes are presented in Figure 6 and Tables 10 and 11. The improved accuracy of the new scheme is observed.

**Figure 6.** Comparison of different moments with selection function *<sup>S</sup>* = *xy*, and breakage functions *b*3, *b*4.

**Table 10.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = *xy* and breakage function *<sup>b</sup>*<sup>3</sup> at *<sup>t</sup>* <sup>=</sup> 10.



**Table 11.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* = *xy* and breakage function *<sup>b</sup>*<sup>4</sup> at *<sup>t</sup>* <sup>=</sup> 10.

#### *4.2. Three-Dimensional Examples*

**Volume conservative problems:** In this section, we consider two test problems with size-dependent selection function *<sup>S</sup>*(**y**) <sup>=</sup> *<sup>φ</sup>*(**y**). The three-dimensional extension of the above-mentioned breakage functions *b*<sup>1</sup> and *b*2, that is,

$$b\_{\mathbb{S}}(\mathbf{x}|\mathbf{y}) = \frac{2}{\psi(\mathbf{y})'} \quad \text{and} \quad b\_{\delta}(\mathbf{x}|\mathbf{y}) = 2\delta\left(\mathbf{x}\_1 - \frac{y\_1}{2}\right)\delta\left(\mathbf{x}\_2 - \frac{y\_2}{2}\right)\delta\left(\mathbf{x}\_3 - \frac{y\_3}{2}\right)$$

are considered here. Like before, we can only calculate the zeroth and first moments exactly and they are given in the following Table 12.

**Table 12.** Exact moment functions for *<sup>S</sup>* <sup>=</sup> *<sup>φ</sup>*(**y**) and breakage functions *<sup>b</sup>*5, *<sup>b</sup>*6.


From Figure 7 and Tables 13 and 14, we can observe that the WFV scheme estimates the moment functions more accurately as compared to the EFV scheme.

(**a**) Breakage function *b*<sup>5</sup>

**Figure 7.** *Cont*.

**Figure 7.** Comparison of zeroth- and first-order moments with selection function *<sup>S</sup>* = *<sup>φ</sup>*(**x**), and different breakage functions *b*5, *b*6.

**Table 13.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* <sup>=</sup> *<sup>φ</sup>*(**x**) and breakage function *<sup>b</sup>*<sup>5</sup> at *<sup>t</sup>* <sup>=</sup> 10.


**Table 14.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* <sup>=</sup> *<sup>φ</sup>*(**x**) and breakage function *<sup>b</sup>*<sup>6</sup> at *<sup>t</sup>* <sup>=</sup> 10.


**Hypervolume conservation:** In this instance, we consider the size-independent daughter distribution function *<sup>b</sup>*7(**x**|**y**) <sup>=</sup> <sup>8</sup> *<sup>ψ</sup>*(**y**) along with the constant selection *<sup>S</sup>* <sup>=</sup> 1 and size-dependent selection *<sup>S</sup>*(**y**) = *<sup>ψ</sup>*(**y**). The exact moments are calculated in the Table 15.

**Table 15.** Exact moment functions for *<sup>S</sup>*(**y**) <sup>=</sup> 1, *<sup>S</sup>* <sup>=</sup> *<sup>ψ</sup>*(**y**) and breakage function *<sup>b</sup>*7.


Figure 8 and Table 16 exhibit the improved accuracy obtained from the WFV scheme over the standard scheme to predict the above mentioned three moments.

**Figure 8.** Comparison of zeroth- and first-order moments with selection function *<sup>S</sup>* <sup>=</sup> 1, and breakage functions *b*7.

**Table 16.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* <sup>=</sup> 1 and breakage functions *<sup>b</sup>*<sup>7</sup> at *<sup>t</sup>* <sup>=</sup> 10.


On the other hand, Figure 9 and Table 17 represent the comparison of the above moments as predicted by the WFV and EFV schemes.

**Figure 9.** Comparison of zeroth and cross moments with selection function *<sup>S</sup>* = *<sup>ψ</sup>*(**x**), and breakage function *b*7.

**Table 17.** Relative error for the weighted moments at different grid points for the test case with *<sup>S</sup>* <sup>=</sup> *<sup>ψ</sup>*(**x**) and breakage function *<sup>b</sup>*<sup>7</sup> at *<sup>t</sup>* <sup>=</sup> 10.


#### **5. Conclusions**

In this article, we have proposed finite volume schemes for solving multidimensional fragmentation problems. In addition to the one-dimensional scheme of [23], it is also extended in the multidimensional setup. It is observed that a careful reconstruction of the standard multidimensional scheme leads to the development of very accurate schemes. The newly proposed schemes obey the conservation laws and also predict several physical moment functions with high accuracy. Several empirical test problems in two- and threedimensions are collected from the literature to validate the efficiency of the new schemes.

**Author Contributions:** Conceptualization, J.S. and A.B.; methodology, J.S.; software, J.S.; validation, J.S. and A.B.; formal analysis, J.S.; investigation, J.S.; resources, J.S and A.B.; data curation, J.S.; writing—original draft preparation, J.S.; writing—review and editing, J.S. and A.B.; project administration, J.S. and A.B.; funding acquisition, J.S. and A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** J.S. thanks the NITT seed grant (NITT/R&C/SEED GRANT/19-20/P-13/MATHS/JS/E1) for their funding support during this work. A.B. gratefully acknowledges funding by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation)—Project-ID 416229255—SFB 1411. The authors acknowledge support by Deutsche Forschungsgemeinschaft and Friedrich-Alexander-Universität Erlangen-Nürnberg (FAU) within the funding program Open Access Publishing.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Acknowledgments:** The authors are thankful to the anonymous reviewers for providing valuable suggestions which have helped to improve the quality of the manuscript.

**Conflicts of Interest:** The authors declare no conflict of interest.

### **Appendix A. Proof of Conservation Laws in Two-Dimensions**

We first calculate the difference <sup>F</sup>*i*+1/2 − F*i*−1/2.

$$\begin{split} \mathcal{F}\_{i+1/2,j} - \mathcal{F}\_{i-1/2,j} &= \sum\_{l\_1=i+1}^{I} \sum\_{m\_1=j}^{J} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1}^{\mathcal{S}} \sum\_{l\_2=1}^{i} \mathcal{B}\_{l\_2,j[l\_1,m\_1} \Delta\_{l\_2} \Delta\_{l\_1,m\_1} \\ &- \sum\_{l\_1=i}^{I} \sum\_{m\_1=j}^{J} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1}^{\mathcal{S}} \sum\_{l\_2=1}^{i-1} \mathcal{B}\_{l\_2,j[l\_1,m\_1} \Delta\_{l\_2} \Delta\_{l\_1,m\_1} \\ &= \sum\_{l\_1=i}^{I} \sum\_{m\_1=j}^{J} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1}^{\mathcal{S}} \mathcal{B}\_{i,j[l\_1,m\_1} \Delta\_{i} \Delta\_{l\_1,m\_1} \\ &- \sum\_{m\_1=j}^{I} n\_{l,m\_1} \mathcal{A}\_{i,m\_1}^{\mathcal{S}} \sum\_{l\_2=1}^{i} \mathcal{B}\_{l\_2,j[l\_2,m\_1} \Delta\_{l\_2} \Delta\_{i,m\_1} . \end{split} \tag{A1}$$

In a similar manner, calculating and simplifying the difference <sup>G</sup>*i*+1/2 − G*i*−1/2.

$$\begin{split} \mathcal{G}\_{i,j+1/2} - \mathcal{G}\_{i,j-1/2} &= \sum\_{l\_1=i}^{I} \sum\_{m\_1=j+1}^{J} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1}^a \sum\_{m\_2=1}^{j} \mathcal{B}\_{i,m\_2|l\_1,m\_1} \Delta\_{m\_2} \Delta\_{l\_1,m\_1} \\ &- \sum\_{l\_1=i}^{I} \sum\_{m\_1=j}^{J} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1}^a \sum\_{m\_2=1}^{j-1} \mathcal{B}\_{i,m\_2|l\_1,m\_1} \Delta\_{m\_2} \Delta\_{l\_1,m\_1} \\ &= \sum\_{l\_1=i}^{I} \sum\_{m\_1=j}^{J} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1}^a \mathcal{B}\_{i,j|l\_1,m\_1} \Delta\_{j} \Delta\_{l\_1,m\_1} \\ &- \sum\_{l\_1=i}^{I} n\_{l\_1,j} \mathcal{A}\_{l\_1,j}^a \sum\_{m\_2=1}^{j} \mathcal{B}\_{i,m\_2|l\_1,j} \Delta\_{m\_2} \Delta\_{l\_1,j} . \end{split} \tag{A2}$$

Next, we calculate the following flux,

H*i*+1/2,*j*+1/2 − H*i*+1/2,*j*−1/2 − H*i*−1/2,*j*+1/2 − H*i*−1/2,*j*−1/2 = *I* ∑ *<sup>l</sup>*1=*i*+<sup>1</sup> *J* ∑*<sup>m</sup>*1=*j*+<sup>1</sup> *nl*1,*m*1A*l*1,*m*<sup>1</sup> *i* ∑ *<sup>l</sup>*2=<sup>1</sup> *j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*l*1,*m*1Δ*l*2,*m*2Δ*l*1,*m*<sup>1</sup> − *I* ∑ *<sup>l</sup>*1=*i*+<sup>1</sup> *J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*l*1,*m*<sup>1</sup> *i* ∑ *<sup>l</sup>*2=<sup>1</sup> *j*−1 ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*l*1,*m*1Δ*l*2,*m*2Δ*l*1,*m*<sup>1</sup> − *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> J* ∑*<sup>m</sup>*1=*j*+<sup>1</sup> *nl*1,*m*1A*l*1,*m*<sup>1</sup> *i*−1 ∑ *<sup>l</sup>*2=<sup>1</sup> *j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*l*1,*m*1Δ*l*2,*m*2Δ*l*1,*m*<sup>1</sup> + *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*l*1,*m*<sup>1</sup> *i*−1 ∑ *<sup>l</sup>*2=<sup>1</sup> *j*−1 ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*l*1,*m*1Δ*l*2,*m*2Δ*l*1,*m*<sup>1</sup> . = <sup>−</sup> *I* ∑ *<sup>l</sup>*1=*i*+<sup>1</sup> *nl*1,*j*A*l*1,*<sup>j</sup> i* ∑ *<sup>l</sup>*2=<sup>1</sup> *j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*l*1,*j*Δ*l*2,*m*2Δ*l*1,*<sup>j</sup>* + *I* ∑ *<sup>l</sup>*1=*i*+<sup>1</sup> *J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*l*1,*m*<sup>1</sup> *i* ∑ *<sup>l</sup>*2=<sup>1</sup> B*l*2,*j*|*l*1,*m*1Δ*l*2,*j*Δ*l*1,*m*<sup>1</sup> + *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> nl*1,*j*A*l*1,*<sup>j</sup> i*−1 ∑ *<sup>l</sup>*2=<sup>1</sup> *j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*l*1,*j*Δ*l*2,*m*2Δ*l*1,*<sup>j</sup>* − *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*l*1,*m*<sup>1</sup> *i*−1 ∑ *<sup>l</sup>*2=<sup>1</sup> B*l*2,*j*|*l*1,*m*1Δ*l*2,*j*Δ*l*1,*m*<sup>1</sup> . (A3)

Further rearrangement and simplification of the terms gives

$$\begin{split} \mathcal{H}\_{i+1/2,j+1/2} - \mathcal{H}\_{i+1/2,j-1/2} - \mathcal{H}\_{i-1/2,j+1/2} - \mathcal{H}\_{i-1/2,j-1/2} \\ = \sum\_{l\_1=1}^{I} \sum\_{m\_1=j}^{I} n\_{l\_1,m\_1} \mathcal{A}\_{l\_1,m\_1} \mathcal{B}\_{i,j[l\_1,m\_1} \Delta\_{i,j} \Delta\_{l\_1,m\_1} - \sum\_{l\_1=i}^{I} n\_{l\_1,j} \mathcal{A}\_{l\_1,j} \sum\_{m\_2=1}^{j} \mathcal{B}\_{i,m\_2[l\_1,j]} \Delta\_{i,m\_2} \Delta\_{l\_1,j} \\ - \sum\_{m\_1=j}^{I} n\_{l\_1,m\_1} \mathcal{A}\_{i,m\_1} \sum\_{l\_2=1}^{i} \mathcal{B}\_{l\_2,j[i,m\_1} \Delta\_{l\_2,j} \Delta\_{i,m\_1} + n\_{i,j} \mathcal{A}\_{i,j} \sum\_{l\_2=1}^{i} \sum\_{m\_2=1}^{j} \mathcal{B}\_{l\_2,m\_2[l\_1,j]} \Delta\_{l\_2,m\_2} \Delta\_{i,j} . \end{split} \tag{A4}$$

Therefore, substituting relations (A1), (A2) and (A4) in the discrete formulation (34) and simplifying, we get

d*ni*,*j*Δ*i*,*<sup>j</sup>* d*t* = *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*<sup>β</sup> l*1,*m*<sup>1</sup> B*i*,*j*|*l*1,*m*1Δ*i*,*j*Δ*l*1,*m*<sup>1</sup> − *J* ∑ *<sup>m</sup>*1=*<sup>j</sup> ni*,*m*1A*<sup>β</sup> i*,*m*<sup>1</sup> *i* ∑ *<sup>l</sup>*2=<sup>1</sup> B*l*2,*j*|*l*1,*m*1Δ*l*2,*j*Δ*l*1,*m*<sup>1</sup> + *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*<sup>α</sup> l*1,*m*<sup>1</sup> B*i*,*j*|*l*1,*m*1Δ*i*,*j*Δ*l*1,*m*<sup>1</sup> − *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> nl*1,*j*A*<sup>α</sup> l*1,*j j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*i*,*m*2|*l*1,*m*1Δ*i*,*m*2Δ*l*1,*m*<sup>1</sup> − *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> J* ∑ *<sup>m</sup>*1=*<sup>j</sup> nl*1,*m*1A*l*1,*m*1B*i*,*j*|*l*1,*m*1Δ*i*,*j*Δ*l*1,*m*<sup>1</sup> <sup>+</sup> *I* ∑ *<sup>l</sup>*1=*<sup>i</sup> nl*1,*j*A*l*1,*<sup>j</sup> j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*i*,*m*2|*l*1,*j*Δ*i*,*m*2Δ*l*1,*<sup>j</sup>* + *J* ∑ *<sup>m</sup>*1=*<sup>j</sup> ni*,*m*1A*i*,*m*<sup>1</sup> *i* ∑ *<sup>l</sup>*2=<sup>1</sup> B*l*2,*j*|*i*,*m*1Δ*l*2,*j*Δ*i*,*m*<sup>1</sup> − *ni*,*j*A*i*,*<sup>j</sup> i* ∑ *<sup>l</sup>*2=<sup>1</sup> *j* ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*i*,*j*Δ*l*2,*m*2Δ*i*,*<sup>j</sup>* = *I* ∑ *<sup>l</sup>*1=*i*+<sup>1</sup> *J* ∑*<sup>m</sup>*1=*j*+<sup>1</sup> *nl*1,*m*1A*l*1,*m*1B*i*,*j*|*l*1,*m*1Δ*i*,*j*Δ*l*1,*m*<sup>1</sup> − *ni*,*j*A*i*,*<sup>j</sup> i*−1 ∑ *<sup>l</sup>*2=<sup>1</sup> *j*−1 ∑ *<sup>m</sup>*2=<sup>1</sup> B*l*2,*m*2|*i*,*j*Δ*l*2,*m*2Δ*i*,*j*. (A5)

Temporal evolution of total volume: Taking sum over *i* and *j* of Equation (A5), we get

$$\begin{split} \frac{\mathbf{d}}{\mathbf{d}t} [\mathcal{M}\_{1,0} + \mathcal{M}\_{0,1}] &= \sum\_{i=1}^{I} \sum\_{j=1}^{J} \sum\_{l\_{1}=i+1}^{I} \sum\_{m\_{1}=j+1}^{J} n\_{l\_{1},m\_{1}} \mathcal{A}\_{l\_{1},m\_{1}} \mathcal{B}\_{i,j|l\_{1},m\_{1}} \Delta\_{i,j} \Delta\_{l\_{1},m\_{1}} \\ &- \sum\_{i=1}^{I} \sum\_{j=1}^{J} n\_{i,j} \mathcal{A}\_{i,j} \sum\_{l\_{2}=1}^{i-1} \sum\_{m\_{2}=1}^{j-1} \mathcal{B}\_{l\_{2},m\_{2}|i,j} \Delta\_{l\_{2},m\_{2}} \Delta\_{i,j} .\end{split}$$

Changing the order of summation in the first term, we get

$$\begin{split} \frac{d}{dt}[\mathcal{M}\_{1,0} + \mathcal{M}\_{0,1}] &= \sum\_{i=1}^{I} \sum\_{l\_{1}=i+1}^{I} \sum\_{m\_{1}=1}^{I} \sum\_{i=1}^{l\_{1}-1} \sum\_{j=1}^{m\_{1}-1} n\_{l\_{1},m\_{1}} \mathcal{A}\_{l\_{1},m\_{1}} \mathcal{B}\_{i,j|l\_{1},m\_{1}} \Delta\_{i,j} \Delta\_{l\_{1},m\_{1}} \\ &- \sum\_{i=1}^{I} \sum\_{j=1}^{J} n\_{l,j} \mathcal{A}\_{i,j} \sum\_{l\_{2}=1}^{i-1} \sum\_{m\_{2}=1}^{j-1} \mathcal{B}\_{l\_{2},m\_{2}|i,j} \Delta\_{l\_{2},m\_{2}} \Delta\_{i,j} \\ &= 0. \end{split}$$

Hence, the volume conservation law is obeyed.

Temporal evolution of zeroth moment: Dividing both sides of Equation (A5) by *xi* <sup>+</sup> *yj* , and taking sum over *i*, *j*, we get

$$\begin{split} \frac{\mathbf{d}\mathcal{M}\_{0,0}}{\mathbf{d}t} &= \sum\_{i=1}^{I} \sum\_{j=1}^{I} \frac{1}{(\mathbf{x}\_{i} + \mathbf{y}\_{j})} \sum\_{l\_{1}=i+1}^{I} \sum\_{m\_{1}=j+1}^{I} n\_{l\_{1},m\_{1}} \mathcal{A}\_{l\_{1},m\_{1}} \mathcal{B}\_{i,j|l\_{1},m\_{1}} \Delta\_{i,j} \Delta\_{l\_{1},m\_{1}} \\ &- \sum\_{i=1}^{I} \sum\_{j=1}^{J} \frac{n\_{i,j}\mathcal{A}\_{i,j}}{(\mathbf{x}\_{i} + \mathbf{y}\_{j})} \sum\_{l\_{2}=1}^{i-1} \sum\_{m\_{2}=1}^{j-1} \mathcal{B}\_{l\_{2},m\_{2}|i,j} \Delta\_{l\_{2},m\_{2}} \Delta\_{i,j} .\end{split}$$

Changing order of summation in the first term, we get

$$\begin{split} \frac{\mathbf{d}\mathcal{M}\_{0,0}}{\mathbf{d}t} &= \sum\_{l\_{1}=1}^{I} \sum\_{m\_{1}=1}^{J} n\_{l\_{1},m\_{1}} \mathcal{A}\_{l\_{1},m\_{1}} \Delta\_{l\_{1},m\_{1}} \sum\_{i=1}^{l\_{1}-1} \sum\_{j=1}^{m\_{1}-1} \left[ \frac{1}{\mathbf{x}\_{i} + \mathbf{y}\_{j}} - \frac{1}{\mathbf{x}\_{l\_{1}} + \mathbf{y}\_{m\_{1}}} \right] \mathcal{B}\_{i,j|l\_{1},m\_{1}} \Delta\_{i,j} \\ &= \sum\_{l\_{1}=1}^{I} \sum\_{m\_{1}=1}^{J} \frac{n\_{l\_{1},m\_{1}} \mathcal{A}\_{l\_{1},m\_{1}}}{\Phi\left(\mathbf{x}\_{l\_{1}}, \mathbf{y}\_{m\_{1}}\right)} \Delta\_{l\_{1},m\_{1}} \sum\_{i=1}^{l\_{1}-1} \sum\_{j=1}^{m\_{1}-1} \left[ \frac{\Phi\left(\mathbf{x}\_{l\_{1}}, \mathbf{y}\_{m\_{1}}\right)}{\Phi\left(\mathbf{x}\_{i}, \mathbf{y}\_{j}\right)} - 1 \right] \mathcal{B}\_{i,j|l\_{1},m\_{1}} \Delta\_{i,j} .\end{split}$$

#### **Appendix B. Exact Moments**

The temporal evolution of zeroth moment <sup>M</sup>0,0(*t*) is given by

$$\frac{\mathrm{d}\mathcal{M}\_{0,0}(t)}{\mathrm{d}t} = \int\_{0}^{\infty} \int\_{0}^{\infty} \frac{1}{\mathrm{x} + y} \left[ \frac{\partial \mathcal{F}(t, \mathbf{x}, y)}{\partial \mathbf{x}} + \frac{\partial \mathcal{G}(t, \mathbf{x}, y)}{\partial y} - \frac{\partial^{2} \mathcal{H}(t, \mathbf{x}, y)}{\partial \mathbf{x} \partial y} \right] \mathrm{d}y d\mathbf{x}.\tag{A6}$$

For *<sup>S</sup>*(*x*, *<sup>y</sup>*) = 1 and *<sup>b</sup> x*, *y*|*x* , *y* = 2 *<sup>x</sup> <sup>y</sup>* , we have

$$\frac{\partial \mathcal{F}(t, \mathbf{x}, \mathbf{y})}{\partial \mathbf{x}} = -\int\_0^\mathbf{x} \int\_y^\infty \frac{\varepsilon + y}{\mathbf{x} + \upsilon} \frac{2}{\mathbf{x}\upsilon} n(t, \mathbf{x}, \upsilon) \mathrm{d}\upsilon \mathrm{d}\varepsilon + \int\_\mathbf{x}^\infty \int\_y^\infty \frac{\mathbf{x} + \mathbf{y}}{\mathbf{u} + \upsilon} \frac{2}{\mathbf{u}\upsilon} n(t, \mathbf{u}, \upsilon) \mathrm{d}\upsilon \mathrm{d}\mathbf{u} \mathrm{d}\upsilon$$

$$\frac{\partial \mathcal{G}(t, \mathbf{x}, \boldsymbol{y})}{\partial \boldsymbol{y}} = -\int\_{\mathbf{x}}^{\infty} \int\_{0}^{y} \frac{\boldsymbol{x} + \boldsymbol{\xi}}{\boldsymbol{u} + \boldsymbol{y}} \frac{2}{\boldsymbol{u} \boldsymbol{y}} n(t, \boldsymbol{u}, \boldsymbol{y}) \mathrm{d}\boldsymbol{\xi} \mathrm{d}\boldsymbol{u} + \int\_{\mathbf{x}}^{\infty} \int\_{\boldsymbol{y}}^{\infty} \frac{\boldsymbol{x} + \boldsymbol{y}}{\boldsymbol{u} + \boldsymbol{v}} \frac{2}{\boldsymbol{u} \boldsymbol{v}} n(t, \boldsymbol{u}, \boldsymbol{v}) \mathrm{d}\boldsymbol{v} \mathrm{d}\boldsymbol{u} \mathrm{d}\boldsymbol{v}$$

and

$$\begin{split} \frac{\partial^2 \mathcal{H}(t, \mathbf{x}, y)}{\partial \mathbf{x} \partial y} &= \int\_0^\chi \int\_0^y \frac{\mathbf{c} + \xi}{\mathbf{x} + y} \frac{2}{\mathbf{x} \eta} n(t, \mathbf{x}, y) \mathrm{d}\mathbf{x} \mathrm{d}\mathbf{c} - \int\_\mathbf{x}^\infty \int\_0^y \frac{\mathbf{x} + \xi}{\mathbf{u} + y} \frac{2}{\mathbf{u} \eta} n(t, \mathbf{u}, y) \mathrm{d}\xi \mathrm{d}\mathbf{u} \\ &- \int\_0^\chi \int\_y^\infty \frac{\mathbf{c} + \chi}{\mathbf{x} + \upsilon} \frac{2}{\mathbf{x} \upsilon} n(t, \mathbf{x}, \upsilon) \mathrm{d}\upsilon \mathrm{d}\varepsilon + \int\_\mathbf{x}^\infty \int\_y^\infty \frac{\mathbf{x} + \chi}{\mathbf{u} + \upsilon} \frac{2}{\mathbf{u} \upsilon} n(t, \mathbf{u}, \upsilon) \mathrm{d}\upsilon \mathrm{d}\mu. \end{split}$$

Now substituting in Equation (A6) and simplifying, we get

$$\begin{split} \frac{\mathbf{d}\mathcal{M}\_{0,0}(t)}{\mathbf{d}t} &= \int\_{0}^{\infty} \int\_{0}^{\infty} \frac{1}{\mathbf{x}+\mathbf{y}} \left[ \int\_{x}^{\infty} \int\_{y}^{\infty} \frac{2(\mathbf{x}+\mathbf{y})}{\mathbf{u}\mathbf{v}(\mathbf{u}+\mathbf{v})} n(t,\boldsymbol{\mu},\boldsymbol{\nu}) \mathbf{d}\mathbf{v} d\mathbf{u} \right. \\ &\quad \quad \quad \quad \quad \quad \quad \quad - \int\_{0}^{\mathbf{x}} \int\_{0}^{y} \frac{2(\mathbf{e}+\mathbf{y})}{\mathbf{x}\mathbf{y}(\mathbf{x}+\mathbf{y})} n(t,\mathbf{x},\boldsymbol{\mu}) \mathbf{d}\mathbf{f} d\mathbf{e} \right] \mathbf{d}\mathbf{y} d\mathbf{x} \\ &= \int\_{0}^{\infty} \int\_{0}^{\infty} \int\_{0}^{u} \int\_{0}^{v} \frac{2}{\mathbf{u}\mathbf{U}} \mathbf{d}\mathbf{x} d\mathbf{y} f(t,\boldsymbol{\mu},\boldsymbol{\nu}) \mathbf{d}\mathbf{v} d\mathbf{u} - \int\_{0}^{\infty} \int\_{0}^{\infty} f(t,\mathbf{x},\boldsymbol{y}) \mathbf{d}\mathbf{y} d\mathbf{x} \\ &= \mathcal{M}\_{0,0}(t). \end{split}$$

Solving these equations, we get

$$
\mathcal{M}\_{0,0}(t) = \mathcal{M}\_{0,0}(0) \exp(t). \tag{A7}
$$

Similarly, one can calculate the exact moment functions corresponding to the different selection and breakage functions.

#### **References**

