POLYDISPERSE GRANULAR FLOWS OVER INCLINED CHANNELS

This thesis focuses on modelling the dynamics of dense granular materials flowing over an inclined channel, utilising a continuum description. In the process of understanding and developing this, besides continuum modelling, the thesis also exhibits the utility of discrete particle simulations (DPMs), and the need for developing an accurate micro-macro mapping technique.

As most of these dense gravity-driven flows are shallow in nature, for monodisperse mixtures, Chapter 2 illustrates the formulation of a novel one -- dimensional (width- and depth-averaged) shallow granular model. Using this model, we not only predict the flow dynamics -- flow height, velocity and granular jumps or shocks -- but also shows that one can forecast the existence of multiple steady states for a given a range of channel inclinations. However, in reality, the majority of flowing particulate mixtures are known to comprise of particles with varied physical attributes, i.e. they are bidisperse or polydisperse. Thereby, as a step towards understanding the associated flow dynamics, and developing improved continuum models, several studies presented in this thesis have utilised discrete particle method. DPMs provide a plethora of information at a particle scale, such as particle position, velocity, interaction forces or stresses. In order to accurately map the particle scale mechanics onto a macroscopic continuum scale, Chapter 3 comprehensively presents a generic framework for an efficient and accurate micro-macro mapping technique for polydisperse mixtures comprising of convex shaped particles, e.g. spheres. More importantly, the method presented is valid for any discrete data, e.g. particle simulations, molecular dynamics and experimental data, for both steady and unsteady configurations.

Before employing the efficient mapping technique of Chapter 3 to its full capacity, based on the current understanding of bidisperse segregation dynamics, we formulate in Chapter 4 a mixture theory segregation model for bidisperse mixtures varying both in size and density. The developed formulation is an extension to an already existing size-segregation model, and is applicable to both shallow (linear velocity profile) and thick (Bagnold profile) flows. Besides predicting the extent of segregation, the theory also predicts zero or weak segregation for a range of size and density ratios, which was further benchmarked using DPMs. Although, we developed an efficient continuum size- and density-segregation model, a detailed study is to be implemented in order to determine more accurate pressure scalings and further extend it to polydisperse mixtures. As a stepping stone, towards determining these pressure scalings, in Chapter 5 we give an example application of the micro-macro mapping technique (illustrated in Chapter 3). For simplicity, we consider a purely size-based segregation model, which was built upon a pressure scaling function containing an unknown parameter. Not only did we determine this unknown material parameter but, more importantly, we also found out that the complete size- and density-based segregation in any flowing particulate mixture is an effect of the generated kinetic stress, rather than the contact stress. The current form of the existing scaling functions is, however, still an active area of research, which definitely needs further attention and care.

Chapters 3, 4 and 5, show how one can mix and match continuum models with DPMs using an efficient coarse-graining method. However, it is still vital to see if the DPMs can actually emulate reality. As a consequence, we illustrate in Chapter 6, how DPMs can be used as a suitable alternative to experiments using two commonly used DPM experiments.


Granular materials
Everything seems simpler from a distance.

-Gail Tsukiyama
Matter exists in a variety of forms, of which granular materials -although simplebehave differently from any of the other standard forms (solid, liquid and gas). For example, a sand castle appears like a solid, which when tumbled upon, crumbles down like an avalanche. On a few occasions, this avalanche takes down the whole castle, whereas sometimes it is confined only to a thin layer on its surface (liquid). In another example, shake up some crushed ice in a cocktail mixer and it behaves like a gas, whereas when trying to pour some salt out of an orifice of a jar, the flow tends to choke and gets clogged at the orifice. Granular material has it all, solid, liquid, gas, plastic flow, glassy, etc. -it can mimic all of these behaviours. Thus exhibiting a prime reason for investing our interest in understanding these simple but, mysterious materials.

History
Our knowledge of granular matter dates back to an era even before mankind's existence. These materials exist on a range of scales -micron sized to planetary -in a variety of forms, exhibiting multitudes of interesting phenomena. The earliest encounter between man and granular matter happened in the surroundings in which he evolved, i.e. with soil, sand, stones, rocks, etc., but more importantly it was his need for survival that ultimately led to an inevitable bond between humanity and granular matter, existing till today. Historical records indicate that primitive men often settled along the river banks where fresh food was available throughout the year. However, due to climatic changes coupled with migration of primitive tribes, new food gathering methods had to be devised. Man needed to find some way to store the nutrients for periods of time when no fresh food was available. Seeds, such as cereal grains, seemed to be one way to solve Introduction the storage problems. Ultimately around 8000 BCE, during neolithic revolution, this seed storage led man towards becoming more agriculturally oriented. Grains were likely among the first cultivated crops. They could be grown relatively easily, farmed in surplus quantities and were suitable for storage in harsh cold winter climates. Thereby, marking the dawn of grain agriculture -man's preliminary experience with granular materials.
Since the neolithic revolution, man began to understand the importance of granular materials. Besides agricultural advances, newer innovations, such as wheeled wagons (to transport), pottery (to store), ploughs (to cultivate), etc., came into light. The economy flourished due to trade and manufacturing of these food resources (grains). By 6000 BCE the emergence and spread of food productions established the social and economic foundations for a civilisation. With the innovation of sun-baked bricks from clay, sand and water mixtures, civilisations arose in the period 3800 BCE -3500 BCE. New cities were built acting as trade, economic, political and religious centres. Granular materials began to play a crucial role in the economic development; the scale and nature of the trade had now changed, with substantial amounts of precious materials, such as lapis lazuli, and huge quantities of everyday commodities, such as grain, textiles, timber and metal ores, being imported to and exported from major cities. As time evolved, human ability to manipulate granular materials, such as sand, stones, etc., advanced qualitatively. Not only did it contribute towards urbanisation of civilisations, but also led to the construction of mega structures, such as the pyramids of Giza (circa 1700 BCE), the lighthouse of Alexandria (circa 280 BCE) and many more. And since then, knowingly or unknowingly, granular materials have become an integral component of our daily life.
As years progressed, new empires arose and perished with borders being replotted every few years. While mankind was and still is busy fighting over territorial or racial or 1 Introduction between micro-metre to centimetre scale, solely due to easy manufacturing, and availability of experimental measurement techniques. The majority of this thesis focuses on mixture particles of different sizes and densities. Furthermore, interactions between individual particles are dissipative in nature, i.e. energy or momentum exchange between them involve inelastic mechanisms. When in contact, e.g. due to a collision, energy is lost in several forms such as heat, sound, plastic deformation, etc. For example, if a box of plastic balls is shaken or vibrated, the energetic balls eventually come to rest, thus proving the dissipative nature of these materials. In this thesis, we consider rigid, dry grains and do not consider soft particles, cohesive effects, or interaction with a surrounding fluid, making dissipation as one of the major factors that distinguishes granular matter from fluids.
In nature and many industrial settings -during flows in rotating drum mixers, over conveyor belts, on inclined channels feeding materials into hoppers, silos, etc. -granular media often experience a variety of external forces, such as collisions, shaking, shear, compaction, etc. When subjected to these forces, not only do the dynamics -mixture state, velocity, forces, etc. -of these materials evolve, but they also yield a galore of interesting phenomena, such as particle segregation, granular Leidenfrost and many more states. As an example, let us consider a model dry granular mixture containing only two differently sized marbles. When allowed to flow over a sufficiently long rough inclined channel, besides the varying flow dynamics -flow height and velocity, granular jumps or shocks -the larger marbles end up near the surface whereas the smaller ones settle near the base of the flow. Similarly, if the same mixture is rotated a few times in a thin horizontal drum mixer, we observe the smaller marbles moving radially inwards to form a central core, while the larger marbles move radially outward, surrounding the core. This phenomenon is defined as particle segregation due to differences in size, which, for example, the majority of industries would like to prevent, as an inhomogeneous mixture blend hampers their product quality. In reality, however, granular mixtures often comprise of particles with differences in several of their physical attributes rather than just size alone. Thereby, making them more complicated, challenging and most importantly an interesting area of research.
Since decades, granular media has been a subject of many studies, ranging from static conditions to flowing, dry hard-particles to soft-particles immersed in a fluid. This thesis focusses on dynamical systems, such as the rapid flow of dense 2 dry spheres over a rough inclined channel. We make an attempt to understand and predict the motion or the dynamics of such gravity-driven granular flows through continuum models, discrete particle simulations and, more importantly, by utilising accurate discrete to continuum mapping methods, i.e. the micro-macro transition. Before we progress with the descriptions of our findings and results, in the following chapters, we proceed by briefly introducing the methods we employ for modelling these rapid dense granular flows.

Modelling dry dense granular flows
In the years before 1980, the majority of studies [2] -concerning granular mixturesfocussed on industrial settings. Although extensive, most of the problem solving procedures relied on small scale pilot projects, which often looked for instant solutions, and came up with a list of do's and dont's, e.g. see Johanson [3], based on empirical findings. There definitely existed very little theoretical work concerning these flows. Thus costing invaluable time and money. As years progressed, technological advances allowed for a much improved and systematic comprehension of flowing granular materials over inclined channels. Based on this understanding, modelling of these dense granular flows has seen some remarkable developments, both on theoretical and numerical fronts.

Continuum theories
Owing to several industrial (e.g. mining, pharmaceutical, food processing, etc.) and geophysical applications (e.g. landslides, pyroclastic flows, debris flow, etc.), development of continuum formulations to model cohesionless dense granular flows down an inclined channel came into light in the late 80's. These are briefed in the following subsections.

Exploiting the shallowness
Often on daily news channels, we come across hazardous events that constitute a big threat to human safety. Some examples include the Elm rockfall (Switzerland) in 1881 [4][5][6], the Sherman Glacier rock avalanche (Alaska) [7], the prehistoric Blackhawk landslide (Alaska) [8] and the various ice avalanches that keep disassociating themselves into motion. All these landslides, rockfalls and snow and ice avalanches that dislodge themselves off the steep slopes, travel large distances before they come to rest. Given the risks these geophysical events pose at human safety and property, predicting them would not only help one avoid catastrophic damages, but also allow for designing efficient safety measures.
One of the characteristic feature of these flows is that they are shallow in nature, i.e. the ratio of the characteristic length scales associated with the flow depth H to the ones associated with the downstream flow length L is small, H /L << 1. Using this shallowness argument, studies formulated continuum models -describing flows over inclined channels -which are classified into two categories; shallow-water-like theories extended to granular flows (hydraulic avalanche models) by Grigorian et al. [9] (and others [10][11][12]) and the Mohr-Coulomb type model put forward by Savage and Hutter [13], which have further been extended in [14][15][16][17][18] to account for both complex basal topography and interstitial pore fluid. Thus, allowing one to mathematically describe the granular avalanches.
The Savage-Hutter model can be systematically derived from the general mass and momentum balance laws. The model uses a simple Coulomb's sliding friction law at the base and the model is closed by assuming that the granular material is always in a stress state consistent with the Mohr-Coulomb yield criterion. On exploiting the shallowness of the flow, the leading-order mass and momentum balance laws are integrated through the flow depth to obtain a one-dimensional theory along the flow direction for the avalanche thickness and the downslope velocity. A depth-averaged valuef of a vari-Introduction able f is given byf = 1 h h 0 f d z. Thus, the resulting equations in Cartesian coordinates, dropping the bars, are then where d is the rate of deposition, θ is the local angle of inclination of the slope (basal topography), h(x, t ) and u(x, t ) is the local flow depth and downslope velocity. The source term (driving force) S is given by S = (tan θ − µ) cos θ, (1.2) where µ is the Coulomb sliding friction coefficient. The Earth-pressure coefficient, K , arises from the Mohr-Coulomb yield criterion. The shape factor, α 1 =ū 2 /(ū) 2 , arises from the depth-averaging. This model has been utilised to quantitatively predict the spreading of granular material flowing over inclined channels. At first glance, the Savage-Hutter theory has a strikingly analogous structure when compared to the shallow-water equations [19]. However, the constitutive properties significantly complicate the model by introducing a highly nonlinear Earth-pressure coefficient K into the theory, which pre-multiplies the pressure in the downslope momentum balance. For K = 1, the Savage-Hutter model is reduced to the shallow-water-like model of Grigorian et al. [9]. Thus assuming that the flowing granular material acts like a shallow inviscid fluid with a Coulomb friction law. Furthermore, both models assume the granular material to be incompressible. This thesis employs the Savage-Hutter model. For more details regarding the derivation of the shallow granular model, see Bokhove and Thornton [20]. Thence, with suitable closure laws determined, the above shallow granular model is able to quantify the possible flow quantities associated with any cohesionless granular material flowing over an inclined channel. Many geophysical studies, however, report of complex flow deposits or patterns that arise due to the multi-component -differences in particle size, density, shape etc. -aspect of granular materials. These differences have a predominant effect on the bulk dynamics. Geophysical findings report that when a granular material avalanches down the underneath rough topography, an inverse grading in size is observed, i.e. the small sized particles settle near the bottom and the larger sized particles rise towards the surface. When this inversely graded material is further sheared, large particles tend to migrate towards the front of the flow and smaller ones towards the rear. This can have significant effects on the bulk dynamics as larger particles are susceptible to greater resistance than the smaller ones. Due to these phenomena, the bulk flow may tend towards forming lobate fingers [21,22] or spontaneously self-channelise to form lateral levees that enhance the run-out distance of, e.g., debris flows [23].
Thus, to account for these additional aspects, suitable accurately predicting segregation models need to be formulated, which when combined with the above shallow granular model would allow one to quantify these complex phenomena. For example, Woodhouse et al. [24] developed a model for these finger formations by coupling 1 7 the above shallow granular model with a depth-averaged size-based particle segregation model. Furthermore, the friction coefficient in ( 1.2) was considered to be a function of the local volume fraction of the small particles. However, the particle segregation model used is applicable to bidisperse mixtures varying in size alone. Thus illustrating the need to continue develop more segregation models which can take into account polydisperse mixtures, varying particle densities and shapes, etc. This brings us to the following section, where we briefly describe the existing segregation models.

Segregation models
The British Materials Handling Board [25] holds particle segregation to be responsible for non-uniform mixture blends, which in turn lead to poor product quality and processing difficulties.
Studies have proposed several mechanisms [26] to be responsible for particle segregation due to differences in their physical properties, such as size [27], density [28], shape [29], inelasticity [30], surface roughness and friction [31]. However, differences in size and density are the primary factors for de-mixing in dense free-surface flows over inclined channels. Although, buoyancy effects appears to explain the reasons for density-based segregation, there still exist a debate concerning the mechanism responsible for size-based segregation. For years, kinetic sieving [32,33] and squeeze expulsion [34] were regarded as the dominant mechanism for size-based segregation. In this mechanism, Savage and Lun [34] proposed the concept of a random fluctuating sieve. The basic idea behind kinetic sieving is that, as the granular material avalanches downslope of an inclined channel, smaller-sized particles have a larger probability to fall into the gaps that open up beneath them. This is complemented by squeeze expulsion, which levers all the particles in the upwards. Thus resulting in a net flux of small particles towards the base and large grains towards the free-surface of the flow. Recently, Fan and Hill [35] proposed an alternative hypothesis for dense systems, where large particles are driven to regions of higher velocity fluctuations by kinetic stress gradients, which results in them rising to the surface of the avalanche in a rotating drum [36] or being driven to the sidewalls in a vertical chute flow [35]. The problem can also be viewed as one of lift and drag forces acting on a large particle in an effective fluid medium of fine particles [28,37]. This thesis, however, focuses on kinetic sieving and squeeze expulsion. In addition, the thesis focusses on the effects of both size and density differences. Although, density-based segregation is weaker than kinetic sieving, it is still strong enough to prevent particle-size segregation altogether [38], if the large particles are sufficiently dense, and promotes size-segregation when the small particles are denser.
Based on this understanding of percolation and diffusion, given x, y and z are the downslope, cross-slope and depth direction, Bridgwater et al. [39] were the first to formulate a continuum model quantifying particle segregation in a bidisperse mixture of particles varying in size alone. Their equation governing the granular mixture state is in terms of the volume concentration, φ, of a component expressed as a fraction of the solid volume, ∂φ ∂t where t is the time, z is the downslope direction, q and D are the percolation and diffu-Introduction sion rate. However, they soon realised that the rate of percolation, q, was dependent on the shear rate, the particle size ratio and the normal pressure. As years progressed, Savage and Lun [34] used statistical mechanics and information entropy theory to arrive at a segregation model from first principles. Their model was formulated in terms of number densities and fluxes. Although the model from Savage and Lun [34] considered various functional forms for the shear rate, it certainly had a downside because the model predicted segregation even in the absence of gravity, which is odd given kinetic sieving is a gravity driven process. From a different perspective, Dolgunin and Ukolov [40] developed a model on the basis of an equivalent mass transfer equation, which accounts for the granular mass transfer due to convection, quasi-diffusion and segregation, ( 1.4) where u is the downslope velocity. Although, the above model ( 1.4) has all the features essential for describing particle segregation, a general framework to derive such models was still lacking. In year 2005, Gray and Thornton [41] proposed this general framework by utilising the principles of mixture theory [42] (as also utilised in this thesis). The theory states that each constituent can simultaneously occupy both space and time, resulting in overlapping partial fields, and satisfies the mass and momentum balance laws, which are stated in terms of these partial fields as below ∂ρ ν ∂t + ∇ · ρ ν u ν = 0, ∂ ∂t ρ ν u ν + ∇ · ρ ν u ν ⊗ u ν = ∇ · σ ν + ρ ν g + β ν .
( 1.5) The superscript 'ν' denotes the constituent type. Variables ρ ν , u ν , σ ν denote the partial density, velocity and stress corresponding to each constituent type-ν. The variable g represents the gravity vector and β ν is the force experienced by each constituent due to other constituents. The partial stress tensor σ ν is considered to be a sum of two components, the volumetric −p ν 1 and the deviatoric τ ν stress, such that σ ν = −p ν 1 + τ ν . By using the shallowness argument, it can be shown that the pressure dominates in the normal direction, and both the deviatoric stresses and normal acceleration terms can be neglected. They further observed that as the small particles percolate downwards they support a smaller part of the load when compared to the larger particles. Using these assumptions after some manipulations and scalings, including a suitable form for the drag force β ν , see 4, they [41] arrived at a mixture theory segregation model, where u, v, w is the downslope, cross-slope and normal (depth-direction) velocity, respectively, S r is a dimensionless segregation rate and F (φ) = φ(1 − φ) is the driving segregation flux. The above model was further extended to consider particle diffusion [43], polydisperse mixtures [44] and higher order flux functions [45,46]. However, the above stated models considered size-based segregation alone, except [45] see Chap. 4. Recently, by introducing the particle size as an independent coordinate, Marks et al. [47] Introduction have led to innovative optimisations that have helped reduce the high computational time. Recent studies have been focussing on a GPU 3 -based framework for developing a highly parallelised GPU based DPM solver, see for a brief performance overview for CPU-and GPU-based systems, see Chapter 6 and references therein. With such massive parallelism now available using programmable GPU hardware, rapid advances have been witnessed during the past 10 years. Innovative simulations on the order of millions of particles are being conducted in quasi real-time in confined environments [60], rotating drums [61], blenders [62], granular soils [63]. This is because of the sophisticated algorithms developed for efficient collision detection [64][65][66][67][68][69][70][71] and constructing memory-efficient data structures [72]. Using high performance GPU's, [73] investigated the size effects in granular mixture flows, whereas [74] used it to simulate fractures in heterogeneous media. Additionally, to bridge the gap between ideal and realistic mixtures, studies have also considered to simulate non-spherical particles using the GPU-based framework, see [75] for triangular particles, [76,77] for convex polyhedrals.
On the whole, despite the need of calibrating and validating DPMs, particle simulations have shown to be an efficient, much appreciated, alternative to experiments.

Micro-macro transition
Given that DPM is an efficient tool to be utilised to probe the intricate details of granular dynamics, there still lies a gap between the discrete and continuum models. Macroscopic continuum quantities, such as density, velocity, stresses and other necessary fields are essential in any analysis involving validation or calibration of, e.g., a continuum model.
The mapping of the microscopic scale dynamics onto a macroscopic continuum scale has been under focus since the classical studies by [78,79] and others [80]. Based on a variety of theoretical postulates, various methods for micro-macro transition have been formulated to extract these macroscopic quantities efficiently, e.g. binning of the microscopic fields into small volumes [81] and the method of planes [82]. However, most of them are restricted in terms of their application due to various limitations, see [81] and the references therein. One of the challenges or requirements for multi-scale methods is to efficiently map the microscopic particle dynamics onto a macroscopic fields, which in turn satisfy the classical equations of continuum mechanics, i.e. the fundamental balance law of mass and momentum Dρ + ρ∇ · u = 0, D(ρ u) + ρ u∇ · u = ∇ · σ + ρ g . (1.7) The above equations are stated in terms of the mass density ρ, bulk velocity u, and stress tensor σ. Coarse graining approaches to granular materials first appeared in the work of [83] and has been extended by various studies [84][85][86][87][88][89][90][91][92][93][94][95]. The coarse graining techniques have two essential advantages over other types of averaging techniques. These are (i) the macroscopic quantities exactly satisfy the continuum laws of motion and (ii) they are applicable to both static and dynamic granular media. With these advantages, the coarse graining approach has been utilised to study the results of computations or experiments 3 Graphics processing unit 1 11 and their characterisation in terms of density, velocity, stress, strain, couple-stress and other fields; for 2D granular systems [96][97][98][99][100], hopper flows [101,102]. Furthermore, the coarse graining method described in [95] has been extended to granular mixture flows near boundaries or discontinuities [103,104] and bidisperse mixtures [105]. These extensions [103] have been been applied to analyse shallow granular flows [56] and segregation phenomena in bidisperse granular mixtures [105].

Thesis outline
The fundamental goal of this thesis is to utilise and extend the aforementioned theories.
As a stepping stone, in Chapter 2, we begin by considering shallow monodisperse granular mixture flowing over a rough inclined channel. In addition, the inclined channel also has contracting sidewalls located downstream near the channel exit. On exploiting the shallowness argument in the cross-slope direction, we further width-average the depthaveraged shallow granular equations. Thus, resulting in a novel one-dimensional granular hydraulic theory. For simplicity, we assume the Earth-pressure coefficient K = 1 and the shape factor α 1 = 1. In addition, the simple Coulomb-like friction law is replaced by a much more efficient empirically determined constitutive law, validated using discrete particle simulations. Given this, flow profiles are predicted in a very simple and efficient manner. As a verification step, the solution of the one-dimensional model is compared with an equivalent two-dimensional shallow granular model. However, no validation with discrete particle simulations or experiments is carried out. With a goal towards analysing bidisperse mixture flows using both continuum theories and discrete particle simulations, Chapter 3 focusses on extending the efficient micro-macro mapping technique to multi-component mixtures. Using the concepts of mixture theory, coarse graining expressions for macroscopic partial quantities are systematically constructed. In an attempt to test the limits of these coarse graining expressions, we apply them on discrete particle data obtained from the simulation of a bidisperse mixture (varying in size alone) for both steady and unsteady scenarios. The derived expressions are generic and can be easily extended to multi-component mixtures varying in both size and density.
Chapter 4 showcases a classic example of using discrete particle simulations as a tool to validate a developed continuum model. Using the same principles of mixture theory, an existing bidisperse purely size-based segregation model is extended to take into account the density differences as well. However, no diffusive remixing is taken into account. By doing so, the resulting theory predicts zero segregation for a range of size and density ratios. As a logical step and as an alternative to intensive experiments, we utilise the discrete particle simulations (no micro-macro mapping), to validate the theoretical prediction.
Given we have the coarse graining (CG) expressions at hand, Chapter 5 illustrates a simple but effective application of these expressions. In a simple mixture theory sizeand density-based segregation model, unknown parameters arise based on the flux functions used (1.6). Using the particle data set as utilised in Chapter 4 and the CG expressions from Chapter 3, macroscopic continuum fields are constructed. Using these fields closure parameter is determined as required to complete the continuum model.
Finally, in Chapter 6 we present an extensive review where we advocate discrete par-

One-dimensional shallow granular model
We consider monodisperse dry granular mixture flowing down an inclined channel, with a localised contraction: theoretically and numerically, using the shallow granular theory. For closure, we consider an empirically determined and discrete particle simulation validated constitutive friction law, which also accounts for the existence of steady uniform flows for a range of channel inclinations. From the depth-averaged shallow granular theory, we present a novel extended one-dimensional granular hydraulic theory, which for steady flows predicts multiple flow regimes like smooth flows without jumps or steady shocks upstream of the channel or in the contraction. For supercritical flows, the onedimensional model is further verified by solving the two-dimensional shallow granular equations through a discontinuous Galerkin finite element method. On comparison, the one-and two-dimensional solution profiles, averaged across the channel width, surprisingly match although the two-dimensional oblique granular jumps largely vary across the converging channel.
This chapter is still in prep. for a publication.

Introduction
A considerable number of industrial processes involve materials in a granular form. Grains of dissimilar properties are often mixed, fed, or separated through a variety of devices in these processes. Partially filled rotating drums and blenders are used in pharmaceutical and food production industries [1], whereas rotary kilns and inclined cylinders [2] are associated with chemical processes involving sinter 1 , cement and iron production due to their ability for continuous material feed. In Europe, [3] sales of granular material, which at some stage of a production process is poured, mixed, or separated, is estimated to be over 6 billion euros with a production of over one million tonnes annually. Amongst several particle transport mechanisms associated with industrial processes, this work concerns analysing rapid free surface granular flow similar to what is observed in steel manufacturing at any global steel company. At a steel manufacturing site, the iron production process involves the inflow of sinter, pellets and coke, via a rotating inclined channel, into a blast furnace for iron-ore melting; see Fig. 2.1. Here, pellets are spheres of certain diameter produced by shaping finely ground particles of iron ore. The mixture is pelletised as it is easy to feed, produce and store. As the sinter, pellet and coke mixture is fed in a layered pattern into the blast furnace; non-uniformity in the mixture properties like size and density leads to particle segregation. Furthermore, the rough base of the hopper is fitted with rivets, which makes the base uneven and rougher. Thus, leading to complex flow and deposition patterns and, more importantly, implying that a thorough understanding of the dynamics of such complex flow phenomena is essential in improving the iron quality, control of the production process (efficiency) and the design of the devices handling mixture feed. As a stepping stone, this chapter considers the flow of a monodisperse mixture alone, i.e. mixture comprised of same type of particles, which also implies no particle segregation, over non-rotating rough inclined channels with constrictions.
In reality, the majority of granular flows in nature (avalanches, landslides, etc.) and industries dealing with inclined channel flows are shallow, i.e. the ratio of the characteristic length scales in the normal (H ) to the streamwise direction (L) is small, H /L << 1. Although qualitative understanding of monodisperse mixture flows over inclined channels has existed for some time; several avalanche models, by exploiting the shallowness aspect, proved to be successful in quantitatively analysing these granular flows. In essence, an avalanche model utilises the already existing shallow water theory from the fluids community and extends it to model shallow granular free-surface flows. However, one needs to know the corresponding granular constitutive relations required to relate the normal and the tangential stresses. To our knowledge, the earliest known extension of the shallow water theory was implemented by Grigorian et al. [4], to predict the snow avalanche paths in the Ural mountains. The formal existence of shallow granular (SG) theory was established by Savage and Hutter [5], who averaged the mass and momentum balance equations in the depth direction (depth-averaging) and assumed a Mohr-Coulomb rheology with a constant Coulomb basal friction law. In depth-averaging, one averages out the depth dependency from the flow quantities, such as the flow height and velocity. For more details also see Bokhove and Thornton [6]. With this established set of shallow granular equations, Savage and Hutter [7] and Greve and Hutter [8] further extended the SG theory, to predict the flow of an initially stationary finite mass of cohesionless granular material down a variable basal topography, both concave and convex. As years progressed, studies extended and generalised SG theory to consider flows over varying topography, see [9][10][11][12][13][14], and used them for predicting precarious zones in the alps (a mountain range in Europe), e.g., see [15][16][17]. Further simplifications were made by Gray et al. [18], where the in-plane deviatoric stresses is assumed to be negligible. This assumption reduces the set of equations of Savage and Hutter [5] to bear a superficial resemblance to the nonlinear shallow water equations with source terms. This reduced set of hyperbolic equations has been successfully used to precisely predict granular flows past obstacles (e.g. pyramids [18], wedges [19] ), through constrictions [20,21] and cylinders [22]. Besides theoretical applications, studies have also utilised them to design avalanche deflecting walls for the protection purpose, e.g. see [21,23]. We focus here on utilising the shallow granular theory for inclined channel flows through contractions.
In Vreman et al. [20], alongside listing several other approaches attempted to predict and understand granular dynamics; the motion of granular matter flowing over a smooth inclined channel, constrained by contracting sidewalls, was investigated by means of theoretical, numerical and experimental analysis. Results revealed upstream moving bores or shocks, a stable reservoir state and weak oblique shocks. Flow states and flow regimes were explained via a unique granular "hydraulic" theory described by a set of equations obtained by extending the asymptotic analysis of Gray et al. [18] to onedimension. Instead of granular material, Akers and Bokhove [24] analysed water flow on a horizontal plane constrained downstream by contracting sidewalls, theoretically and experimentally. In this chapter, the same one-dimensional hydraulic theory is extended to the granular case, including frictional effects. These one-dimensional shallow water equations combined with a theoretical and experimentally determined constitutive friction law leads to a one-dimensional shallow granular continuum model. To close this model we use a friction law (constitutive law/closure relation) determined by Pouliquen and Forterre [25] and verified through discrete particle simulations of Weinhart et al. [26]. The discrete particle simulations construct a map between micro-scale and macroscale variables and functions thereby determining the closure relations needed in the continuum model.
Through the one-dimensional asymptotic model, the flow regimes observed for granular flow in an inclined channel -with a contraction -are illustrated on a F 0 − B c plane, where F 0 is the channel upstream Froude number and B c is the critical nozzle width ratio. The latter is defined as the ratio of the nozzle width W c and the upstream channel width W 0 . The one-dimensional asymptotic theory gives a thorough, although approximate overview of the possible flow regimes. Results obtained via the one-dimensional asymptotic theory are later verified by solving the two-dimensional shallow granular equations using a discontinuous Galerkin finite element method (DGFEM). Solving the two-dimensional shallow granular equations through DGFEM not only helps in verification of the asymptotic theory, but most importantly it scrutinises whether the constitutive friction law holds in higher space dimensions (two-dimensions). Finally, although not shown in this chapter, the results from the continuum model should be compared Shallow granular model with full DPM simulations in order to investigate the validity of the assumptions of the granular shallow-layer model.

Asymptotic theory
In rapid free surface shallow granular flows, the ratio of the characteristic length and velocity scales in the normal to streamwise direction is small (<< 1). By utilising the asymptotic analysis from Gray et al. [18], along with a series of approximations listed in Appendix A, we obtain the depth-averaged shallow granular equations. The dimensional 2D shallow granular model is stated below (2.1) The above 2D shallow granular equations, (2.1), are derived via a step by step procedure presented in the book chapter, Bokhove and Thornton [27]. Eqns. (2.1) represent the conservation of mass and momentum in terms of the flow quantities, i.e. flow depth h = h(x, y, t ) and velocity (u, v) = (u(x, y, t ), v(x, y, t )) in a channel of width W = W (x) with basal topography b(x, y), where x and y is the down-and cross-slope direction, respectively. Variable t denotes time, µ(h, u) is the basal friction coefficient and K a material constant denoting stress anisotropy. The subscripts t , x, and y denote the respective partial derivatives and g is the acceleration due to gravity. The variable θ denotes the chute angle, which is chosen such that the average inter-particle and particle-wall forces are in balance with the downstream force of gravity acting on granular particles, leading to a uniform flow in the absence of a contraction. Using the aspect ratio argument as used for depth-averaging, the flow quantities across the channel are also averaged, see Appendix A. 1. After averaging the flow quantities across the chute, i.e. averaging out the y-dependence, the 2-dimensional model equations (2.1) reduce to the 1-dimensional depth-and width-averaged shallow granular equations The flow quantities h and u are independent (averaged out) of y, i.e. h = h(x, t ), u = u(x, t ), with g n = g cos θ. We consider a uniform channel with a constant upstream width W 0 and a localised linear contraction, which is monotonically decreasing in width, W (x) ≤ W 0 , from x m till it reaches the minimum nozzle width W (x c ) at x = x c . Here, x m and x c are, respectively, the x location of the contraction entrance and exit, Fig. 2

.2.
We always consider the minimum nozzle width at the end of the channel. Thence, we introduce the dimensionless variables denoted by primes, The variables W 0 , u l and h l are suitable characteristic length scales for the flow domain (usually a reference value for the channel width), flow velocity and flow depth, respectively. In our model, W 0 , u l and h l are the values defined upstream of the contraction at x = −x l < (x m = 0). We also define an upstream Froude number as F l = u l / g n h l . The average slope of the contraction is given by, α = (W 0 − W c )/(x c − x 0 ). After substituting the scaled variables and dropping the primes, we obtain a non-dimensional depth-and Shallow granular model width-averaged shallow-layer model for isotropic (simpler case), K = 1, granular flow, For a detailed derivation, see Appendix A. 1. The 1-dimensional shallow-layer equations (2.4) constitute of the continuity equation and the downslope momentum equation. One could also arrive at these equations via a standard control volume analysis of a column of granular material viewed as a continuum from the base to the free surface, using Reynolds-stress averaging and a leading order closure. Moreover, in order to have a closed system of shallow-layer equations we need a constitutive friction law, which we shall briefly describe in the following section.

Constitutive law/Closure relation
The basic difference between the shallow-layer fluid model and a granular one is the presence of the basal friction coefficient µ, where µ is the ratio of the shear to normal traction at the base. Some of the previously developed dry granular models incorporated a dry Coulomb-like friction law [5]. However, the Coulomb-like friction law holds only in two cases: 1. When the inclined channel is smooth, fully developed uniform flows are found to exist at one critical inclination angle [28][29][30]. Above this angle the material accelerates and below this angle the flowing material eventually stops. The rheological properties of flows over smooth channels are well described by a constant friction constant, which equals the tangent of the angle of friction between the material and the base δ, i.e. µ = tan δ.

2.
Similarly, experimental studies also show that the constant friction coefficient holds for accelerating flows over rough channels at higher inclinations [28,31]. Experimental measurements of the shear forces at the bed show the friction coefficient to be independent of the velocity.
For an intermediate range of angles where steady uniform flows reside [32][33][34], the simple Coulomb friction, however, fails to describe the flow rheology on channels with rough beds. Using accurate experimental measurement techniques Pouliquen [32], Forterre and Pouliquen [35] empirically determined a scaling which allows one to predict the variation in the mean (depth-averaged) velocity as a function of the channel inclination, flow depth and channel roughness, where β and γ are constants and F is the Froude number. The effects of changing the channel roughness, channel inclination and other features like mixture particle size is captured in h st op (θ) without any experimental velocity measurements. The Variable h st op (θ) denotes the critical thickness where the flow arrests or comes to a halt. Each 25 channel inclination has a unique critical thickness, which depends on the channel roughness and particle size, see [32] for more details concerning the measurement of h st op (θ). With this scaling law at hand, Pouliquen and Forterre [25] further expressed the stoppage height, as a function of the angle of inclination, where d is the grain diameter and A is a characteristic dimensionless length scale over which the friction varies. In addition, the above empirical friction law (2.6) is characterised by two angles: the angle at which the material comes to rest δ 1 , below which friction dominates over gravity and the angle δ 2 , above which the material accelerates as gravity dominates friction. It is between these two angles that steady flow is possible.
One can obtain the constitutive friction law on combining (2.5) and (2.6). Using the steady state flow assumption µ = tan θ to hold(approximately) in the dynamic case as well, one obtains an improved, valid for lower Froude numbers, empirical friction law As δ 1 → δ 2 , the Coulomb's model is recovered, see Grigorian et al. [4].

Steady state solutions
By utilising the improved macro-scale constitutive friction law (2.7), the steady flow states in a granular flow -through a contraction -are predicted through the shallowlayer granular model (2.4). We begin by defining a non-dimensional Froude variable as F l takes F l = F 0 for values u 0 , W 0 and h 0 at x = x 0 near the sluice gate or F l = F m for values u m , W 0 and h m at the contraction entrance x = x m . For flows in steady state, from the mass balance equation (huW ) x = 0, we obtain a constant mass flux huW = constant. We define the integration constant as Q and for our scaling we consider Q = 1. The momentum balance equation, in a conservative form, is stated as Using  Eqn. 2.11 is analogous to the equation (7) in [24] An analytical solution is found for a special case where W = W (x) and assuming that for a given steady flow, µ = tan θ (inviscid flow) holds (approximately). Eq. (2.11) can for the inviscid case, be written as which when analytically integrated with respect to x, from the channel upstream position x l to some point channel downstream yields For the given closed system (2.11), the well known critical nozzle condition from Houghton and Kasahara [36], F = 1, plays the role of a boundary condition at the channel exit. At this condition the flow at the nozzle is "sonic" or "critical", such that the flow speed u equals the gravity wave speed h/F l , which in terms of dimensional quantities means that the flow speed u equals g n h. On utilising this critical condition for the inviscid case, the solution with F (x c ) = 1 and for x < x m , where F l = F 0 , is 14) The average solutions obtained are smooth as long as the flow is subcritical with F 0 < 1 or supercritical with F 0 > 1. The inviscid solution (2.14) divides the F 0 − B c parameter plane into regions characterising the smooth sub-and supercritical flows, see Fig. 2. 3. In order to obtain the demarcating curves for granular flows with frictional effects, we integrate the ordinary differential equation (ODE) (2.11) using the fourth-order Runge-Kutta scheme, known as RK4, starting from the contraction exit x = x c given the critical nozzle condition F (x c ) = 1, B c and the width W = W (x). The location of the contraction exit, x = x c , is given by the relation x c = L cos θ c , with θ c = sin −1 ((W 0 −W c )/2L), where L is the length of the contraction paddle. The Froude number F l and depth h l is prescribed upstream of the channel at x = x l . Either F l = F 0 upstream of the channel or F l = F m at the contraction entrance. Given the critical nozzle condition is at the contraction exit, where F = 1 by definition, a closer look at Eq. (2.11) indicates,  4. Choosing a circle is convenient, as it has an infinite number of slopes and can be easily fitted at the contraction exit x = x c . Moreover, as the radius R → 0, we return back to the initial sidewall geometry (approximately). Once the new contraction exit, x cnew , is determined, such that LHS=RHS, we arrive at a classic limit problem, stated as The limit (2.16), is solved by using the Taylor expansion series, see Appendix A. 1.5. The finite slope F ′ at the new contraction exit x = x cnew is determined for both the inviscid and viscous cases, see Appendix A. 1. 5. On regularisation, the singularity at the contraction exit due to the critical nozzle condition is resolved such that LHS = RHS. For a given sufficiently large F l > 1 or sufficiently small F l < 1; in order to produce the critical (demarcating) curves -for viscid flows -that allows one to distinguish between the smooth super-and subcritical flows and flows with jump, the ODE (2.11) can be integrated in two ways: (i) Non-regularised approach: To avoid the singularity in Eqn. 2.11, the ODE is integrated starting from the initial contraction exit, x c , with Froude number F = lim ε→0 1 ± ε, respectively. Given the singularity is resolved by determining the new contraction exit and the correct slope dF /d x, as described in Appendix. A. 1.5, a regularised approach can now simply use the critical boundary condition F = 1 and integrates the ODE starting from the new nozzle exit x = x cnew . Via either approach, regularised or non-regularised, the ODE (2.11) is integrated to a point upstream, x = x 0 or x = x m , to find a new estimate for F l , since the scaling F l is unknown beforehand as it is part of the solution. Hence, the solution is found iteratively. One can consider the initial value (an educated guess) for F l = F 0 (B c ) obtained from the solution for the inviscid case, here F 0/m is a function of B c . Once the new estimate for F l is obtained, we proceed iteratively until convergence is reached. Thereby we obtain the necessary demarcation curves on the F 0/m − B c plane. In Fig. 2.5, we show the critical curves obtained using both the regularised and non-regularised approach. As illustrated, the constructed solutions are in good agreement. Utilising the same regularisation approach as above, Fig. 2.6 predicts the existence of different flow regimes by plotting the critical curves which demarcate the F 0/m − B c plane into several regions.

Shock solutions
The closed system of equations, (2.4), is hyperbolic and thus can develop discontinuities in the flow field in finite-time. These discontinuities are nothing but jumps in the depth or velocity of the flow that propagate at a well-defined velocity. The jump conditions are derived as follows. We integrate both the non-dimensional depth-and width-averaged mass and momentum balance equation, from X (t ) − δ to X (t ) + δ and take the limit δ → 0. Both h and u are discontinuous at x = X (t ). For simplicity, we define X − as the limit position on the left side of the jump and X + the limit on the right side of the jump, and the shock speed s = −d X /d t . Utilising these definitions, we begin by integrating the continuity equation. On utilising Leibniz's rule, we have Similarly, using the same above approach, we integrate the momentum equation as below As W + = W − = W , the above equation is simplified to  On utilising Leibniz's rule, hu d x = 0 and On further manipulations, the above equation is restated as below Using the jump condition obtained from the continuity equation (2.17), the above equation (2.18) is restated as The above two jump/shock relations are also known as Rankine-Hugoniot conditions. For upstream moving shocks instead of matching the upstream conditions with the nozzle conditions, we relate the depth h u and u u upstream of the shock to h 1 and u 1 downstream of the shock and also to the depth h c and u c at the nozzle end. Using (2.17) and (2.19) with h + = h u , u + = u u and h − = h 1 , u − = u 1 , the jump conditions are restated as Similarly, in order to obtain the jump conditions for a shock in the contraction, we combine the above jump conditions with the Bernoulli equation and mass conservation in the contraction thus leading to Similarly, below we substitute the above scalings in (2.20), to arrive at a few relations useful for further manipulation of (2.21) 1 , (2.24) Furthermore, by using the critical condition (2.22) in (2.21) 2 , we obtain (2.25) On substituting the above relations (2.24)-(2.25) in (2.21) 1 , we systematically proceed to arrive at (2.26) Shallow granular model Thus, through these above manipulations, we arrive at the jump conditions in terms of the upstream Froude number, shock speed and flow depth ratio ( (2.29) The ODE (2.11) is further integrated upstream from F = F m > 1 at x = x m to find our next estimate of F * l at x = x 0 . Generally, F * l = F l , where F l is the scaling used in Eqn. (2.11). Thereby, one begins with the inviscid result F l = F 0 (B c ). In the inviscid case one can use (2.13) with F l = F 1 at the entrance of the contraction and F m = F 0 > 1 to find  35 For comparison purposes, the demarcating curves, obtained either by regularisation or non-regularisation, are almost identical, see Fig. 2. 5. Nevertheless, the regularisation is necessary to make the numerical analysis mathematically sound, i.e without a singularity at F = 1. The super-(F > 1) and subcritical (F < 1) flow profiles are illustrated via F and h versus x plots in Fig. 2. 7. These profiles are obtained by integrating from a point upstream of the channel into the downstream direction to the contraction exit. For flows with granular jumps (discontinuities in flow quantities), the critical condition at the nozzle exit is F = 1. We commence at the nozzle exit and move into the upstream direction. The jumps in the flow quantities are computed by applying the jump condition, and then finding a point where the downstream and upstream profiles match. Thereby, the flow profiles are efficiently computed from the novel granular hydraulic theory. In the following section we shall validate one of these flow profiles (supercritical), using two-Shallow granular model dimensional solutions obtained through DGFEM.

Verification of 1D theory via 2D DGFEM
The asymptotic theory helps to approximately predict the flow regimes in F 0 − B c plane. In order to verify the 1D results, the two dimensional solutions of the shallow granular equations (SGE) are compared with the former. The additional degree of freedom helps us illustrate the results from the 1D analysis in 2D. Moreover, the verification process checks if the constitutive friction law holds in two dimensions. Solutions are obtained by solving the SGE, along with closure relation µ(h, F ), using the discontinuous Galerkin finite element method (DGFEM). The DGFEM is a blend of high resolution finite element and finite volume methods for solving system of equations, as discussed in Di Pietro and Ern [37]. We adopt a second order space discontinuous Galerkin method for the numerical solution of our SGEs. The weak form of the SGE required for using DGFEM is derived below. For more detailed information regarding the space discontinuous Galerkin method one can refer to [38].

Weak formulation for DGFEM
From Eqn. (2.1) after applying the same scaling as in (2.3), the non-dimensional shallowlayer granular flow system of equations can be written concisely in index notation as 32) are the flux and the source term, respectively. The conserved quantities of U are used to model granular bores/jumps. A closed system of equations is obtained by prescribing the initial conditions and the boundary conditions, such as inflow depth or flow velocity and also specifying the friction law µ(h, f ). For details regarding the space elements and tessellation of the domain of our problem, Ω, see Tassi et al. [38]. The weak form of the given closed set of equations (2.31) is obtained by multiplying it with an arbitrary test function W h ∈ V d h and replacing the exact solution U by its approximation (Ω)|v| K ∈ P p (K )}, with P p (K ) as polynomials of degree p on the element K and L 2 (Ω) is the function space of square integrable functions and d = dim(W h ). Given the domain Ω is discretised into k elements K k , by integrating over the space element K k , we have the weak form for our compact system as Applying the Gauss theorem to eq. (2.34) and summing over all the elements results in (2.35) After defining a tessellation T h of the domain Ω, having a boundary ∂Ω, we define Γ as a set of all the faces in T h . If x is a point on a face S in Γ and if n K is the outward unit normal vector at the boundary ∂K k , then the trace of the function V h on the interior of the element boundary ∂K k is defined as and dΓ is an infinitesimal boundary segment of the element K k . Summing over all the elements we could rewrite the boundary integral in Eq. (2.35) as (αF l i j + βF r i j )(n l K j W l hi + n r K j W r hi ) + (n r K j F r i j + n l K j F l i j )(αW r hi + βW r hi )dΓ, with α + β = 1 and α, β ≥ 0. If we consider the flux to be conservative across the element faces such that S F l i j n l K j W l hi dΓ = − S F r i j n r K j W r hi dΓ, (2.36) and n l K j = −n r K j , we get The numerical fluxF i j (U r ,U l , n K ) replaces the flux (αF l i j +βF r i j ) in (2.37). The numerical flux is a function of the traces U l and U r at the element face, where l and r denote the left and right side of the face, respectively. Finally, the weak formulation using the space-DG method is stated for each element as: With the weak formulation at hand we use suitable numerical fluxes like the HLLC numerical flux or the kinetic numerical flux. The discretised formulation arises when the approximated quantities U h and W h are defined in terms of the suitable basis functions. After which we use the explicit third order total variational diminishing (TVD) Runge-Kutta scheme for time integration. The numerical solution is computed through our open-source dgFEM solver, hpGEM.org.

2D DGFEM solutions
The above weak formulation is utilised to implement, in our open-source code, a DGsolver for the 2D shallow granular equations. A contracting-expanding channel is considered in a Cartesian coordinate system (x, y) ∈ [0, 16] x [0, 1]. The channel converges from x = x m = 6 to x = x c = 11 and diverges from x c = 11 till the end of the channel as seen in the top-most illustration of Fig. 2. 8.

Supercritical Flow
For supercritical flows through a contraction, we take as initial conditions the flow height h = 0.5, downstream velocity u = 0.5, velocity in the y-direction v = 0 and the basal topography b = 0. At the inflow boundary, we specify the flow height h = 1, downstream velocity u = 1 and the cross-slope velocity v = 0. At the outflow, all the variables are extrapolated such that U r = U l . The channel sidewalls are considered to have a solid wall boundary condition, as proposed in Ambati and Bokhove [39], imposed as: where t is the unit tangential vector orthogonal to the normal vector n at the sidewall. For a given nozzle width B c and upstream Froude value F L corresponding to the multiple steady state region in Fig. 2.6, the second (middle) illustration of Fig. 2.8 shows oblique weak (not predominant with huge discontinuities) jumps in the converging region of the channel as predicted through the asymptotic theory.

DGFEM vs. 1D asymptotic theory
The 2D solution obtained from the DG-formulation is averaged in the cross-slope direction and compared with the solution profile constructed using the novel one-dimensional granular theory, case (i) in Fig. 2. 7. From Fig. 2.8, one can conclude that the onedimensional asymptotic theory is in a very close agreement with the two-dimensional solution for supercritical flows.

Summary
As a stepping stone towards analysing the monodisperse granular mixture, flowing over rough inclined channels -with constrictions -we use the extensively utilised depthaveraged shallow granular equations. On exploiting the shallowness aspect in the crossslope direction, which is often associated with the flows considered, we width-averaged the depth-averaged shallow granular equations, and arrived at a novel one-dimensional granular hydraulic theory. For closure, we used an empirically determined constitutive friction law, which was also validated using discrete particle simulations by Weinhart et al. [26]. By using the one-dimensional theory, we can compute the flow profiles for any channel opening and inclination in a quick and efficient manner. Besides this ability, we also illustrated -for steady flows -an existence of multiple flow regimes. As a verification step, for supercritical flows alone, the 1D theory is compared to solutions obtained via numerically solving an equivalent depth-averaged shallow granular equations using DGFEM. On comparison, we found a close agreement between the one-and twodimensional theories.

Future work
A complete verification of the 1D theory has to be carried out, i.e. the one-dimensional granular theory needs to be verified for subcritical flows and also for flows with jumps or shocks. Besides the verification of the theory, a thorough validation needs to be performed. Although attempts have been made to validate the above one-dimensional theory using fully three dimensional discrete particle simulations, efforts were so far in vain. Mainly, because the flow profiles in the contraction did not completely match the 1D predictions. However, the flow profiles for flows over an inclined channel without a contraction do agree with solutions obtained using the 1D theory. Thence, suggesting for further study to understand the underlying dynamics in the contraction.

3
From discrete data to continuum fields The method presented in this chapter is valid for any discrete data, e.g. particle simulations, molecular dynamics, experimental data, etc.; however, for the purpose of illustration we consider data generated from discrete particle simulations of bidisperse granular mixtures flowing over rough inclined channels.

Introduction
To formulate accurate continuum models one constantly needs to calibrate and validate them with the available experimental or numerical data, which are discrete in nature. To implement this mapping in an efficient manner, accurate micro-macro transition methods are required to obtain continuum fields (such as density, momentum, stress, etc.) from discrete data of individual elements (positions, velocities, orientations, interaction forces, etc.). This is the focus of this chapter: How to perform the micro-macro transitional step? Many different techniques have been developed to perform the micro-macro transition, from discrete data, including Irving & Kirkwood's approach [2] or the method of planes [3]; we refer the interested reader to [4,5] and references therein. Here, we use an accurate micro-macro transitional procedure called coarse-graining, as described in [4,[6][7][8][9][10][11][12]. When compared with other simpler methods of performing the micro-macro transitions, the coarse-graining method has the following advantages: (i) the resulting macroscopic fields exactly satisfy the equations of continuum mechanics, even near the boundaries, see [10], (ii) the elements are neither assumed to be spherical or rigid, (iii) the resulting fields are even valid for a single element and a single time step, hence no ensemble-averaging is required, i.e. no averaging over several time steps or stamps. However, the coarse-graining method does assume that (i) each pair of elements has a single contact; i.e. elements are assumed to be convex in shape; (ii) the contact area can be replaced by a single contact point, implying that the overlaps are not too large; (iii) the collisions are enduring (i.e. not instantaneous). Often, micro-macro methods employ ensemble-or bulk-averaging to obtain accurate results; therefore, the methods are only valid for homogeneous, steady situations. The coarse-graining method overcomes these challenges by applying a local smoothing kernel, coarse-graining function, with a well-defined smoothing length, i.e. coarse-graining scale, that automatically generates fields satisfying the continuum equations. As an example, one could consider a Gaussian as a coarse-graining function with its standard deviation as a coarse-graining scale. For more details concerning the choice of the coarse-graining functions, see Sec. 3

.2.4.
The coarse-graining method is very flexible and can be used with discrete data from any source, e.g. molecular dynamics, smoothed particle hydrodynamics, discrete particle simulations, experimental data [13], etc. Previously coarse-graining has been successfully extended to allow its application to bulk flows near the boundaries or discontinuities [10,11] and to analyse shallow granular flows [4]. Here, we systematically extend the method to a multi-component unsteady, non-uniform situations, and demonstrate its application by considering the granular flow of spherical particles (convex-shaped). Recently, the technique of coarse-graining was used to analyse steady bidisperse granular mixtures of spheres varying in size alone [43]. Besides extending the technique to unsteady multi-component mixtures, we apply it -for demonstration purpose -to a bidisperse flow of spherical particles, varying in both size and density, over inclined channels for both steady and unsteady configurations. Giving special focus upon the often neglected topic of how to coarse grain in time for unsteady scenarios?
Granular materials, conglomerates of discrete macroscopic objects, are omnipresent, both in industry and nature. Therefore, understanding the dynamics of granular materials [14][15][16] is crucial for a diverse range of important applications, such as predicting natural geophysical hazards [17] to designing efficient material handling equipments [18][19][20][21][22]. Although, in the past 30 years, extensive studies have been carried out in the field of granular materials, today several open-questions in both static and dynamic granular materials are yet to be answered, e.g. failures in static grain silos, rheology of non-spherical flowing grains and many more. In nature, and often in industry, granular materials are polydisperse (multi-component); comprised of elements varying in size, shape, density and many other physical properties [23]. Therefore, in the past few years, much work has been focused on multi-component systems, both experiments and simulations, in a host of different applications, including granular mixture flows in rotating drums [31,32], over non-rotating or rotating inclined channels [33,34], in vibrated beds [35,36], in statics near jamming [37] and many more. Consequently, new continuum models are being formulated that attempt to model the dynamics, e.g. particle segregation, of these multi-facetted granular constituents in different applications [33,[38][39][40][41][42]. In particle segregation, particles often tend to arrange themselves in distinct patterns due to relative differences in their physical attributes. For example, if a bidisperse (twocomponent) mixture -varying in size alone -flows over an inclined channel, eventually the larger particles end up near the free-surface whereas the smaller particles find themselves to appear near the base of the flow [24].
For granular materials, the discrete particle method (DPM) is a very powerful computational tool that allows for the simulation of individual particles with complex interactions [26], arbitrary shapes [27], in arbitrary geometries, by solving Newton's laws for each particle, see [28,29]. Moreover, complex interactions such as sintering, breaking and cohesional particles can be captured, by an appropriate contact model; however, this method is computationally expensive. Nevertheless, with the continuous increase in computational power it is now possible to simulate mixtures containing a few million particles; but, for 1 mm particles this would represent a flow of approximately 1 litre, which is many orders of magnitude smaller than the real life flows found in industrial or environmental scenarios.
Continuum methods, on the other hand, are able to simulate the volume of real environmental and industrial flows, but need simplifying assumptions that often require effective macroscopic material parameters, closure relations or constitutive laws, etc. In order to correctly apply these continuum models, both the continuum assumptions must be validated and the effective material parameters must be determined for a given application; e.g., the Savage-Hutter model [30] for granular geophysical mass flows requires the effective basal friction for closure [4]. However, these continuum models often make assumptions that need to be validated, and contain new continuum properties that must be determined for given materials. These are the so-called validation and calibration steps, which need to be undertaken either by careful experiments or using well chosen small DPM simulations. Thus, motivating the need for an accurate micro-macro method that can deal with multi-component scenarios.

Outline
To extract the averaged macroscopic fields, the coarse-graining (CG) expressions are systematically derived in Sec. 3.2. As a test case, we apply the available CG expressions to bidisperse granular mixtures flowing over an inclined channel, see Fig. 3. 1. In Sec. 3 Figure 3.1: A snapshot of a bidisperse mixture flowing in a periodic box inclined at 26 • to the horizontal (discrete particle simulation). Colours/shades indicate the base/boundary (yellowish green, F b ), species type-1 and type-2 (blue, F 1 and red, F 2 ). We define the bulk as F 1 ∪ F 2 .
for flows in steady state, we show that there exists a range or plateau of smoothing lengths (coarse-graining scale/width) for which the fields are invariant. Although the technique does not require ensemble-averaging, we nevertheless illustrate spatial coarse-graining (averaging in space alone) to be well complemented by temporal averaging (averaging in time). For bidisperse unsteady flows, Sec. 3.3.4 illustrates how to define both spatial and temporal averaging scale such that resolved scale independent time-dependent fields can be constructed. Finally, Sec. 3.4 summarises and concludes our main findings.

Spatial coarse-graining
The current section comprehensively extends the approach of [4,10] to bidisperse spherical systems, and can be easily extended to polydisperse mixtures without any loss of generality. Traditionally, the coarse-graining formulae were derived from the classical laws of conservation of mass, momentum, energy, etc., [9]. Thereby, leading to expressions for the total density, stress, etc., in terms of the properties of all the particles. Here, we generalise this to mixtures (multi-components); therefore, our starting point will be mixture theory [44], which constructs partial mass, momentum and energy balances for each distinct constituent of a mixture.

Mixture theory
As stated above, the coarse-graining formulae will be formulated using the framework of mixture theory, which is often used to study porous media flow problems (e.g. the flow of gas, oil and water mixtures through a deformable porous matrix) [44], sea ice dynamics [45], snow metamorphism [46], determining the properties of concrete [47], swelling of chemically active saturated clays [48] and many more applications.
Mixture theory deals with partial variables that are defined per unit volume of the mixture rather than intrinsic variables associated with the material, i.e. the values one would measure experimentally. The basic mixture postulate states that every point in the mixture is occupied simultaneously by all constituents. Hence, at each point in space and time, there exist overlapping fields (displacements, velocities, densities) associated with different constituents.
Since each constituent is assumed to exist everywhere, a volume fraction Φ ν is used to represent the percentage of the local volume occupied by constituent ν. Clearly, where n is the number of distinct granular constituents in the mixture and φ a denotes the fraction of volume corresponding to interstitial pore space filled with a passive fluid, e.g. air. However, for convenience, studies, e.g. [41], often consider volume fraction of the constituents per unit granular volume rather than per unit mixture volume. As the volume fraction of granular constituents per unit mixture is the volume fraction of each constituent per unit granular volume is which also sum to unity, For each individual constituent, conservation laws for mass, momentum, energy and angular momentum can all be obtained, but here for simplicity, we only consider mass and momentum balance for bulk constituents. Each bulk 1 constituent satisfies the following fundamental laws of balance for mass and momentum [44], (3.5) The above fundamental laws (3.5) are derived from the classical principles of mass and momentum conservation corresponding to each constituent, see [44] for details. ∂ t = ∂/∂t and ∇ = [∂/∂x, ∂/∂y, ∂/∂z] denote the partial temporal and spatial derivatives, respectively. Symbols ' · ' and '⊗' denote scalar and dyadic product. Furthermore, (i) ρ ν and u ν are the partial density and velocity.
(iii) β ν denotes the partial inter-constituent drag force density (drag) which essentially accounts for the net effect of tractions across the interfaces of different constituents. The inter-constituent drag is analogous to the viscous shear tractions resisting the relative motion of fluid through matrix pores.
(iv) b ν represents the partial body force density, which accounts for all the external body forces (generally due to gravity) acting on each constituent ν.
The variables appearing in the theory are partial not intrinsic, these are defined such that their sum is equal to the total mixture quantity. For example, This makes the bulk quantities easy to calculate, by simply summing over all bulk constituents. Various relations can be shown between the intrinsic (by convention a superscript ' * ' denotes an intrinsic variable) and partial variables. In models based on mixture theory, the relationships for velocity and density are For the case where the stress tensor can be represented by a hydrostatic pressure field, it is common in the application of mixture theory [44] to consider a linear volume fraction scaling for the pressure as well, i.e.

A mixture theory for coarse-graining
Consider a DPM simulation with three different types of particles: (bulk) type-1, (bulk) type-2 and boundary, whose interstitial pore-space to is filled with a zero-density passive fluid, see Fig. 3. 1.
will have a radius a i whose centre of mass is located at r i with mass m i and velocity v i . The total force f i (3.9), acting on a particle i ∈ F is computed by summing the forces f i j due to interactions with the particles of the same type j ∈ F ν and other type, j ∈ F /F ν , and body forces b i , e.g., gravitational forces (m i g).
where the Greek subscript α = [x, y, z] denotes the vector components. For each constituent pair, i and j , we define a contact vector Irrespective of the size of constituent i and j , for simplicity, we place the contact point, c i j , in the centre of the contact area formed by an overlap, δ i j , which for small overlaps has a negligible effect on particle dynamics.
To account for the interaction of the two bulk constituents, type-1 and type-2, with the boundary, we will denote the boundary as a third constituent. As the constituents of a bidisperse system are classified under three categories -type-1, type-2, boundary -a three-constituent continuum mixture theory [44] is considered, see Sec. 3.2.1. In other words, we classify the bidisperse system constituents under three categories (i) type-1 constituent (ii) type-2 constituent and (iii) boundary. The set bulk comprising of type-1 and type-2 constituents and F b denotes the boundary constituents, e.g. see Fig. 3.1. Although Fig. 3.1 depicts a flowing (dynamic) system scenario, the above nomenclature is equally applicable to static bidisperse system. For the bulk constituents, F 1 ∪ F 2 , we define partial densities, ρ ν , velocities, u ν , stresses, σ ν , with ν = 1, 2. Additionally, we also define inter-constituent drag force densities, β η→ν , corresponding to the interaction among different constituent with η, ν = 1, 2, b. When η = ν, by definition β η→ν = 0.
(ii) The drag on the bulk constituents due to the boundary is defined as t = β b→1 + β b→2 and is equivalent to the boundary interaction force density (IFD) defined in [10].
In the following sections, using the above postulates of mixture theory, we systematically derive and arrive at the coarse-graining expressions for both partial and bulk quantities in terms of discrete particle data defined above.

Mass density
The partial microscopic (point) mass density for a system (in a zero-density passive fluid) at the point r and time t is given from statistical mechanics as (3.12) where δ( r ) is the Dirac delta function in R 3 . This definition complies with the basic requirement that the integral of the mass density over a volume in space equals the mass of all the particles in this volume.
To extract the partial macroscopic mass density field, ρ ν ( r , t ), the partial microscopic mass density (3.12) is convoluted with a spatial coarse-graining function ψ( r ), see Sec. 3

.2.4, leading to
( 3.13) Essentially, we replace the delta-function with an integrable (real and finite support) coarse-graining function of space, ψ(r ), also known as a smoothing function. For benefits seen later, we define ψ i = ψ( r − r i (t )). From the partial density (3.13), the partial volume fraction is defined as where ρ ν p is the (constant) material density of constituent type-ν. Thereby, the bulk volume fraction is defined as ν = ν 1 + ν 2 . Given the coarse-graining expressions for partial densities (3.13), using (3.11), the bulk macroscopic density field is defined as Thence, on utilising expressions (3.13)-(3.15), one can construct spatially coarse grained fields for partial and bulk density. However, it is still unclear about the choice and type of coarse-graining functions one could use in these expressions. Thereby, in the following section we briefly reflect upon the characteristics and possible forms of coarse-graining functions, ψ( r ).

Which functions can be used to coarse-grain?
The coarse-graining functions ψ( r ) need to possess certain characteristics essential for the technique of coarse-graining: (i) They are non-negative, i.e. ψ( r ) ≥ 0 ensuring the density field to be positive.
(iii) There exists a compact support c ∈ R such that ψ( r ) = 0 for | r | > c.
As a regularisation to the delta-function, below are a selection of archetype cases one could choose from (i) Heaviside: where H represents the Heaviside function and Ω(w) = (4/3)πw 3 is the volume of a sphere in three-dimensional space, with w as its radius.
(ii) Gaussian: fields and is infinitely differentiable. Often a cut-off is utilised in order to compute the fields efficiently.
(iii) Lucy polynomials: In this manuscript, we utilise a family of polynomials called Lucy, see [49]. In three dimensional (3D) space, the 4th-order Lucy polynomial is defined as with c the cut-off radius or the range (compact support) and w = c/2 the coarsegraining scale or pre-determined width (or standard deviation). A Lucy polynomial has at least two continuous derivatives. Moreover, the use of a polynomial form allows one to compute exact spatial averages and gradients of the resulting fields as they are integrable and differentiable analytically.
Note, in all the cases 'w' is defined such that a direct comparison between the different coarse-graining functions for a fixed 'w' can be made.
In the limit w → 0, both the Gaussian and Lucy polynomials tend towards the deltafunction. However, as long as the coarse-graining function is not singular or highly anisotropic, the fields depend only weakly on the choice of the above functions, but strongly on the chosen or predetermined spatial coarse-graining scale, w.
Thus, with the coarse-graining function known and the expressions for partial and bulk mass density at hand, the coarse-graining expressions for partial and bulk momentum density, velocity and stress fields shall be comprehensively derived in the following sections.

Mass balance
By utilising the coarse-graining expression for macroscopic partial mass density (3.13), we derive the governing equation conserving the mass, which is satisfied by each constituent of the mixture. Note that (using the chain rule): is the smoothing kernel around particle i . Using the approach of [9], we consider the time derivative of the coarse-grained partial mass density (3.13). Using (3.17), we have with ν denoting the species-type and p ν ( r , t ) defined as the coarse-grained partial momentum density, The above expression (3.19) corresponds to the microscopic partial momentum density . Moreover, on rearranging the terms in (3.18), using the shorthand notation ∂ t = ∂/∂t and ∇ = [∂/∂x, ∂/∂y, ∂/∂z], we arrive at the mass balance law, in terms of the partial fields, Note that the above result also holds for a single constituent (e.g. single particle) in a mixture, and one does not need to consider an ensemble of constituents, e.g. a collection of particles, to define these fields. Additionally, the macroscopic partial velocity fields, u ν ( r , t ), are defined as the ratios of partial momentum density and mass density fields Thence, the coarse-grained partial mass density and velocity fields are defined such that they exactly satisfy the mixture continuity equation (3.20) which, when summed over the constituent types, leads us to the mass balance law (excluding the boundary) where ρ( r , t ) is the macroscopic bulk mass density field (3.15) and p( r , t ) = ν p ν ( r , t ) is defined as the macroscopic bulk momentum density field. Furthermore, the bulk velocity field, u, is defined as u α = p α ( r , t )/ρ( r , t ), which satisfies the bulk law of mass balance (3.22).

Momentum balance
Besides satisfying mass balance laws, as postulated in mixture theory (Sec. 3.2.1), each constituent (e.g. single particle) of the system also satisfies the fundamental balance law of momentum, which, when stated in terms of partial fields is In order to obtain an expression for the partial macroscopic stress field, σ ν , we rewrite the momentum balance law (3.23) in component form, To begin with, we compute the temporal derivative of p ν α as, is the total force on particle i ∈ F ν . Substituting (3.9), the first term of (3.25) can be expanded as The first term of A ν α , representing interactions between constituents of the same type, satisfies by first interchanging the indices i and j and then applying Newtons' third law, f i j α = − f j iα . On adding the first and the third term from (3.27), it follows that

28) can be restated as
The second term of A ν α , representing inter-species interactions, can be rewritten as

Substituting (3.29) and (3.30) into (3.26), yields
which when simplified results in From the above expression, we define the interspecies drag force density (drag) in (3.24) localised at the contact point c i j . The body force density is defined as To obtain the macroscopic partial stress field σ ν αβ , we use the identity [10] which is rewritten using the chain rule of differentiation and the Leibnitz' rule of integration. In (3.35), b i j = r i − c i j is the branch vector as illustrated in Fig. 3.2. Substituting the expressions (3.35) in A ν α , allows one to compute the force densities along the branch vector between the particles. Using the identity (3.35) and substituting (3.34), A ν α is rewritten as where σ c,ν αβ is the macroscopic partial contact stress field; due to all the contacts among all the constituents. The integral χ i j ensures that the contribution of the force between two constituents i and j to the partial stresses to be proportional to the length of the branch vectors, i.e. the stresses are distributed proportionally based on the fraction of the branch vectors contained within the constituent. Thus, for contacts between a small and a large constituent, the larger sized constituent receives a bigger share of the stress. Following [9], the second term of (3.25), is expressed as where σ k,ν αγ is the macroscopic partial kinetic stress field; Thereby, from (3.39), the total partial stress field, σ ν αβ , is defined as the sum of both partial contact and kinetic stress fields, σ ν = σ c,ν + σ k,ν . Similarly, from (3.10), the total bulk stress field is defined as ( 3.41) In the case of bidisperse mixture, ν = 1, 2, the bulk stress is defined as In order to illustrate a simple application of the above coarse-graining expressions to compute the partial stresses and interspecies drag forces, a simple setup of static bidisperse (large and small) two-dimensional particles (discs) is considered, see Fig. 3. 3. Using the coarse-graining expressions for partial drag (3.34) and stresses (3.39), Fig. 3.3 exhibits the magnitude of partial stresses and drag arising from the contacts between the discs.
So far, we have comprehensively derived and given the coarse-graining expressions for both partial and bulk mass and momentum density, velocity and stress fields including the expressions for the boundary force density, a interspecies drag force density, and the body force density. In the following section, using a convenient medium, we present a simple example to utilise these expressions for a bidisperse mixture where ν = 1, 2. Magnitudes of partial stresses, σ s (small discs type-1) and σ l (large discs type-2), and partial drag experienced by large discs, β l , due to small discs in a static assembly of bidisperse (small and large) twodimensional discs.

Application
Besides the simple example in Fig. 3.3, involving static bidisperse two-dimensional discs, we apply the coarse-graining expressions to a larger bidisperse system in three dimensions (3D). As an example, we consider bidisperse mixtures flowing over inclined channels, as depicted in Fig. 3.1 and described below. This problem was considered previously in [33] and more details of the setup can be found in that article.

Discrete particle simulation (DPM) setup
A fully three-dimensional simulation of an initially homogeneously mixed bidisperse mixture of particles, see Fig. 3.1, is considered. The two different particle types are referred to as type-1 and type-2. If d 1 and d 2 , are defined as the particle diameter of particle type-1 and type-2, then the mean particle diameter is defined as with φ = ν 1 /(ν 1 + ν 2 ) being the volume fraction of particles of type-1.
In our chosen coordinate system, as illustrated in Fig. 3.1, we consider a cuboidal box, set to be periodic in the x-and y-directions and with dimensions (x, y, z) The box is inclined at θ = 26 • and consists of an irregularly arranged fixed particle base, for further details see [4,33] or Chapter 4 or Chapter 5. The parameters in our DPM simulations are non-dimensionalised such that the mean particle diameter d = 1, its mass m = 1 and the magnitude of gravity g = 1 implying the nondimensional time scale t := d /g . The ' ' denotes non-dimensional quantities.
The box is filled with a bidisperse mixture in which the number of particles of each type is where the V box = 20 × 10 × 10 is the volume of the box. The formulae (3.44) ensure that the ratio of total volume of particles of type-1 to the total volume of all the particles is φ and the dimensionless height of the flow, H is the same for all simulations used in this chapter. Using (3.44), for homogeneous initial conditions (randomly mixed), with initial particle volume fraction φ = 0.5, DPM simulations for two different particle size ratios, s = d 2 / d 1 = 2 and 3.5, were carried out.
For the performed simulations, we use a linear spring dashpot model [28,29] with a contact duration of t c = 0.005 d /g , coefficient of restitution r c = 0.88, contact friction coefficient µ c = 0.5 and time step t c / 50. More details about the contact model can be found in [4] and [29].

Spatial coarse-graining
In order to obtain the continuum macroscopic fields, for any stationary or transient particulate system, it is essential to choose a proper spatial coarse-graining scale, w, irrespective of the chosen coarse-graining function, ψ(r ). So the question that arises is how do we choose w? This question is equivalent to asking what do we mean by a continuum description? A continuum description has an implicit length scale associated with it for which the assumptions made in the continuum model are valid and it is this length scale over which we must coarse-grain. When one chooses a length scale, w, smaller than the continuum length scale, the resulting coarse grained data will still show individual particles; these are not continuum fields. On the other hand, if one chooses a large w, it will smear out the macroscopic gradients and the results will be strongly dependent on w. Between these two extremes, their exists a plateau in which the continuum fields obtained are independent of the w chosen and it is this length scale that must be utilised for an efficient micro-macro transition. Thus, leading to another interesting question: Do such plateaus exist for the example we considered?
Quest for the plateaus, i.e. what is an optimal spatial coarse-graining scale?
To determine a suitable scale, bidisperse mixtures of two different particle size ratios s ∈ {2, 3 3.5(b)), respectively. Given these steady flow configurations, we use the above derived coarse-graining expressions to construct the bulk volume fraction, λ(z), as a function of the flow depth, for two different coarse-graining scales, Fig. 3.5(c) ( s = 2.0) and Fig. 3.5(e) ( s = 3.5). By following the steps described in Appendix B, these profiles are constructed by spatially averaging in both x-and y-direction and temporally over a time interval [600, 800] (i.e. 200 snapshots). As seen in these plots, the resulting depth profiles strongly depend upon the chosen coarse-graining scale, w. For s = 2, when averaged on a sub-particle length scale: layering in the flow can be observed near the base of the flow (boundary). However, when averaged on the particle length scale, the layering effect, observed near the base, is smoothened out. The particle-scale density is nearly constant in the bulk, whereas it decays slightly near the base where density oscillations are strong (dilatancy), and near the surface, where the pressure approaches the atmospheric pressure. Thereby, illustrating the larger gradients alone, which are present near the base and the free-surface. The momentum density, velocity and the contact stress show the same qualitative behaviour. Similarly for s = 3.5, The two coloured blocks labelled as '1' and '2' in (d) and (f) denote sub-particle or microscopic scale (1) and particle or continuum scale (2). for a sub-particle length scale, layering is not just observed near the base, but also within the bulk, which is smoothed out when averaged using a particle length scale (denoted by • in Fig. 3.4(f)). However, understanding and illustrating the underlying dynamics of mixtures with larger particle size ratios is beyond the scope of this chapter and will be addressed in a future publication. Nevertheless, an ideal scenario would be to see whether these macroscopic fields are independent of the chosen coarse-graining scale. But, does such a scenario exist? Numerical simulations, see Goldenberg et al. [50] which involve systems of 2D polydisperse disks and Weinhart et al. [51] for monodisperse 3D mixtures flowing over inclined channels, show that for a considerable range of coarse-graining scales, w, the computed fields are independent of the averaging scale. As a step towards our quest for determining this so called range (plateaus), we average these steady state mixture configurations, Fig Thereby, in order to compute meaningful fields for w < 0.01, longer temporal averaging or a larger number of particle ensembles would be needed. In other words implying more particle data needs to be stored, i.e. probably at every 100 (2t c ) time steps. Nevertheless, the existence of this first plateau confirms the presence of a sub-particle length scale, much smaller than the mean particle diameter, for which consistent invariant fields can be defined. We denote this sub-particle scale as microscopic scale. Similarly, for mixtures with particle size ratio s = 3.5, Fig. 3.4(f), the first plateau spans from w = 0.03 -0.2, which is slightly smaller when compared to the one observed in Fig. 3.4(d).
Besides the first plateau, there also exists a second plateau (labelled as 2) in the range of 0.75 ≤ w ≤ 1.5 in Fig. 3.4(d) and 2.3 ≤ w ≤ 3.5 in Fig. 3.4(f). Both plateaus (on particlescale) appear to be narrower than their corresponding first plateaus (effect of using a log-scale for the x-axis). Nevertheless, the presence of the second plateaus confirms the existence of a mean particle length scale for which, again, invariant fields can be constructed. We denote the scales in this range as continuum scale. Moreover, the coarsegraining scales chosen in Therefore, the plots in Fig. 3.4(c-f) show (i) the effects of the chosen spatial coarsegraining scale, w, on the averaging of the fields and (ii) the existence of a range of scales for which invariant fields can be constructed on both sub-particle and particle scale.

Temporal averaging
The choice of a coarse-graining scale for spatial averaging, depends on the scale of the problem, i.e. microscopic or continuum. Now that, for mixtures in steady state, we have determined the ranges/plateaus, from which one could choose a spatial scale, w = w/d , we shift our focus towards investigating the issues concerning temporal averaging of   spatially coarse-grained fields. Thus, leading us to the question: Is spatial averaging complemented by temporal averaging? Note: In the previous section, the fields computed were both spatially and temporally averaged. However, we primarily focussed on the effects of w, the spatial coarse-graining scale, for a fixed temporal averaging width.
In order to carry out in-depth analysis concerning temporal averaging, the same discrete particle simulation as described in Sec. 3. 3.1 is utilised. However, rather than saving data at every 10000 (200t c ) simulation time steps, as done in the previous Sec. 3. 3.2, we consider saving particle data at every 100 (2t c ) simulation time steps, i.e. with the simulation time step d t = 0.0001 (t c /50) we have 100 snapshots for each simulation time unit. For temporal averaging, we consider a fixed averaging time interval, i.e. ∆t a = t min , t max = [652, 1852]. If N a is defined as the number snapshots to average over, for the chosen ∆t a , we have a total of 120000 snapshots. We define this as N a,t ot al .
Given the time interval is defined, we temporally average over N a number of snapshots, which are cleverly chosen from the defined time interval ∆t a ; note that ∆t a = [652, 1852] is fixed. We initially begin with N a = 2 and gradually increase the number of snapshots, N a → N a,t ot al . As a result, for the spatial coarse-graining scale w = 0.1, the effects of N a on temporal averaging of spatially averaged (in x-and y-direction alone) depth profiles of the bulk density are illustrated in Fig. 3.5(b-e). As the value of N a increases, implying an increase in the number of snapshots to average over, the statistical fluctuations gradually disappear, see Fig. 3.5(e). The decrease in these statistical fluctuations due to increasing value of N a can be quantified by computing the L 2 -error, defined as ( 3.45) Note that λ a and λ b are spatially and temporally averaged fields. On plotting E λ against the number of averaging snapshots (N a ), see Fig. 3.5(f), we observe that the error is inversely proportional to the square root of N a , i.e. E λ ∝ 1/ N a , see the dashed line. Finally, from Fig. 3.5, one can infer that, for steady flows, spatial averaging can definitely be complimented by temporal averaging, i.e. there exists an optimal number of snapshots to construct meaningful fields, which in turn is dependent on the chosen spatial coarse-graining scale, w. However, for w > 2.0, effects of the smoothing function take over, leading to overly smooth fields neglecting the boundary effects and their gradients.

Averaging unsteady mixture states
So far, in the previous sections, following the procedure outlined in Appendix B, we have applied our coarse-graining (CG) expressions on particle data corresponding to steady flows 2 . It is, however, the unsteady particle dynamics that is vital for completely understanding the underlying phenomena and developing accurate continuum models. Thereby an essential step would be to examine, in detail, the application of CG expressions to unsteady mixture states.
As an example application, we consider the same system, i.e. of bidisperse granular mixtures (varying in size alone) flowing over inclined channels as described in Sec. 3   For particle size ratio, s = 2, the whole process of segregation happens within the first 500 time units. See Fig. 3.5(a), where the vertical centre of mass, of both large and small particles, is tracked. However, to investigate the application of coarse-graining to transient, unsteady flows, we focus on the part before particle segregation is attained, i.e. when t ∈ [50, 450] see Fig. 3.6(a). Moreover, we consider the dynamics of large particles (partial fields) alone rather than focussing on the bulk. Considering the same data set that was used for our investigation in Sec. 3.3.3 (data stored at every 100 (2t c ) simulation time steps) and following the approach taken in Sec. 3. 3.2, we begin with spatial coarsegraining of particle data available in the time interval ∆t a = [50,450]. As a result, given a spatial coarse-graining scale ( w) is chosen, the spatial averaging is carried out in x-and y-direction alone. Thence resulting in a spatially averaged profile, denoted byζ( t, z In general, given a spatial ( w) and temporal ( w t ) averaging scale, temporal averaging of any spatially averaged field,ζ( t, z), can be defined as where t denotes a point about which we would like to temporally average. Note that: w t determines a time interval over which we temporally average, t − w t , t + w t , see  3.6(d) shows that for a given spatial coarse-graining scale w = 0.4, there exists a range of temporal averaging scales, 30 ≤ w t ≤ 85, for which invariant fields can be constructed. For w t ≥ 90 (N a = 18000), macroscopic averaging (time-smoothening) effects take over and hence leading to a decrease in the density value, whereas for w t < 30, strong statistical fluctuations exist. Similarly, for a given temporal scale, w t = 60 (N a = 12000), the coloured block in Fig. 3.6(e) illustrates that there exists a range of spatial coarse-graining scales for which invariant averaged fields can be constructed, also see w t (data not shown). Thence, implying that there exists a range of both spatial coarsegraining scales and temporal averaging scales for which invariant averaged fields can be computed.
Additionally, we consider a range of spatial w and temporal w t , CG scales, which results in a w t × w phase plot. Thereby, for each combination of a spatial and temporal scale, we spatially and temporally average the available particle data. Once an averaged field is constructed, we track a point, z = 7.0, in the flow depth to analyse its sensitivity to different values of the spatial and temporal scale, similar to what we did earlier. As a result, Fig. 3.7 displays a contour plot for λ L ( z = 7.0) and illustrates that there exists a region of (almost) invariance irrespective of the chosen spatial and temporal averaging scale, see the rectangular region. For w t ≥ 90, macroscopic smoothening effects dominate, while for w t < 30, strong statistical fluctuations exist, as seen in Fig. 3.6(d), and for w > 1.5, effects of large spatial coarse-graining scales take over. Nevertheless, similar regions of invariance are found to be existing at different values of flow depths z and different values of t.
Therefore (i) for a given single data-set, in order to utilise the coarse-graining expressions, see Sec. 3.2, for unsteady flows, one needs to specify both the temporal and spatial scales of averaging, i.e. both spatial and temporal averaging has to be done. (ii) Similar to the results corresponding to steady flows, there exists a range or plateau of temporal and spatial scales for which macroscopic fields can be constructed for flows with timeevolving macroscopic states.

Summary and conclusions
In this work, we comprehensively derived a novel and efficient technique of spatial and temporal mapping, called coarse-graining, for bidisperse systems. The technique can be easily extended to multi-component systems without loss of generality. As an application example, we carried out an in-depth analysis concerning the coarse-graining by using an example bidisperse mixture, of two different size-ratios (same density), flowing over a rough inclined channel, for both steady and unsteady scenarios. Note that this technique is equally applicable to static and polydisperse mixtures as well.
As a result, for steady flows, we have discovered the existence of a range or plateau of spatial coarse-graining scales, both, on the sub-particle (microscopic) and particle (continuum) scale, for which invariant coarse-grained fields can be constructed, see Fig. 3.4. We also found that the spatial averaging is well complemented by temporal averaging, see Fig. 3. 5. Additionally, for unsteady flows, we discovered a region of invariance, see Fig. 3.7, i.e. a range of spatial and temporal coarse-graining scales for which (almost) invariant fields can be constructed.
Here, we did not present any analysis using the coarse-grained quantities to compute the unknown macroscopic parameters [43], or validate continuum formulations and constitutive postulates [4]. This shall be the focus of our future work where we will thrive on developing accurate continuum formulations using the approach of the micromacro transition presented above. Furthermore, no quantitative recommendations are provided as coarse-graining is highly system dependent.
The above coarse-graining method is available as part of an open-source code Mer-curyDPM (mercurydpm.org) and can be run either as a post-processing tool or in real time, see Appendix B. In real-time mode, it not only reduces the data that have to be stored, but also allows for the boundary conditions, etc., to be coupled to the current macroscopic state of the system, e.g. allowing for the creation of pressure-controlled walls.

Mixture theory continuum segregation model
In In the past ten years much work has been undertaken on developing mixture theory continuum models to describe kinetic sieving-driven size segregation. We propose an extension to these models that allows their application to bidisperse flows over inclined channels, with particles varying in density and size. Our model incorporates both a recently proposed explicit formula for how the total pressure is distributed among different species of particles, which is one of the key elements of mixture theory-based kinetic sieving models, and a shear rate-dependent drag. The resulting model is used to predict the range of particle sizes and densities for which the mixture segregates. The prediction of no segregation in the model is confirmed by using discrete particle simulations, and good agreement is found when a single fitting parameter is used which determines whether the pressure scales with the diameter, surface area or volume of the particle.

Segregation model
When free-surface granular flows with particles differing in size and/or density discharge down an inclined plane they often segregate to form complex patterns [2,3]. These flowinduced effects must often be avoided in production processes of the pharmaceutical, chemical, food, iron and cement industry [4,5]. Therefore, a quantitative prediction of segregation is vital in improving the product quality and design of the material handling equipment. Despite its importance, the fundamentals of this phenomenon are incompletely understood.
In general, segregation or de-mixing occurs due to differences in particle properties such as size [6], density [7], shape [8], inelasticity [9], surface roughness and friction [10]. However, differences in size and density are the primary factors for de-mixing in freesurface flows over inclined channels. Experimental studies have been considered to observe the combined effects of size and density difference [11,12] but few continuum models have considered these combined effects [13]. Felix and Thomas [11] experimentally analysed the size and density effects of particles for bi-dispersed mixture flows in rotating tumblers, over inclined channels and pile formations. Using a continuum approach, we present an analysis predicting the degree of segregation in a bidisperse mixture flow, over inclined channels, due to both size and density differences. The zero segregation prediction, from our continuum model, is later benchmarked against the discrete particle simulations.
Among several mechanisms causing segregation [14], our focus lies upon kinetic sieving, which is the dominant mechanism causing segregation in gravity driven free-surface flows, when the size-ratio is less than 2 [15]. For size ratios greater than 2, percolation effects becomes important and should be included in the model [15,16]. We use the framework of mixture theory [e. g. 17] and extend the ideas of Gray and Thornton [18] and Thornton et al. [19] in two ways: (i) we relax the equal particles-species density assumption and (ii) we incorporate a shear rate dependent interspecies drag with a slightly more general pressure scaling function of that proposed by Marks et al. [13]. The resulting theory is able to predict the range of sizes and densities for which segregation will occur directly from the known particle's size and density. Previously Marks et al. [13] stated a way to incorporate density differences, but they have not considered these combined effects in detail.

Particle segregation model.
We choose a domain consisting of a chute inclined at a constant angle θ with respect to the horizontal and a Cartesian coordinate system in which the x-axis points down the chute, the y-axis points across its width and the z-axis points in the upward direction normal to the chute.

Mixture framework.
Our starting point for the model is a granular mixture theory composed of two different constituents, indexed 1 and 2, whose interstitial pore space is filled with air, which has a negligible effect on these dense granular flows. Mixture theory [e.g. 17] for a binary continuum postulates that all constituents of the mixture simultaneously occupy space and time. This leads to overlapping fields with partial pressures, p ν , densities, ρ ν , and 4.1. Particle segregation model.
velocities, u ν = [u ν , v ν , w ν ] T in the three coordinate directions, corresponding to each constituent indexed ν = 1, 2. Each of the constituent satisfies the following fundamental balance laws of mass and momentum for these partial fields where g = (g t , 0, −g n ) T is the gravity vector with g the standard acceleration due to free fall; g t = g sin θ and g n = g cos θ. The β ν 's represent interspecies' drag forces resisting the motion between the constituents. As these forces are internal, from Newtons' third law the sum of these drags must be zero, i.e., β 1 + β 2 = 0. Given a unit mixture volume, each of the constituents occupies a volume fraction φ 1 or φ 2 , including the interstitial pore space. Hence, by definition, the individual volume fractions sum to unity φ 1 + φ 2 = 1. Furthermore, the bulk density ρ, barycentric granular velocity or bulk velocity u and the bulk pressure p are defined as ρ = ρ 1 + ρ 2 , u = (ρ 1 u 1 + ρ 2 u 2 )/ρ , p = p 1 + p 2 respectively. A vital element in the mixture theory is the relation between partial and intrinsic variables. Hereby, variables such as velocity, density and pressure are related as follows with the ' * ' indicating when a variable is intrinsic. Motivated by Marks et al. [13], we assume that f ν scales with the species size s ν as The new definition for partial pressure (4.2) 3 is a slight generalisation of the form used by Marks et al. [13]. When a = 1 or 2 or 3, the pressure precisely scales with the length, surface area, or volume of the particle, respectively.

Drag force.
In a simple shear flow of a mixture with different particle species, the shear causes the species to preferentially fall into the gaps created beneath them (kinetic sieving) or rise upwards through (squeeze expulsion) each other, resulting in interspecies friction [15]. Following Marks et al. [13], a generalised version of Gray and Thornton [18] interaction drag or interspecies' friction is assumed to be with c as a priori unknown coefficient of interspecies' interaction andγ the shear rate. Note that c/γ has a dimension of 1/s. The momentum balance for each individual species (4.1) 2 can be thus restated as In shallow large scale industrial or natural granular flows, the aspect ratio, of velocity and flow length scales in the downslope and cross-slope direction to those in the normal direction, is small. Summing the momentum balance equation (4.5) of each species implies that the flow, at leading order in this aspect ratio, is in lithostatic balance, i.e., ∂p/∂z = −ρg cos θ, see Appendix C. Moreover, at leading order in this aspect ratio, the down-and cross-slope velocity of the species also equals to the bulk down-and crossslope velocity component u ν = u, v ν = v, ν = 1, 2, see Appendix C. Assuming the flow to be in viscous balance by neglecting the inertia terms, substituting the drag force and velocity definitions into the particles' normal momentum balance equation, the species' percolation or normal velocity is (4.6) Here ρ = ρ 2 * /ρ 1 * and s = s 2 /s 1 are the particle density and size ratios, respectively, whereas φ = φ 1 is the volume fraction of species-1. By combining the mass balance equation (4.1) 1 and the percolation velocity (4.6) 1 , we obtain the governing equation for φ 1 = φ as follows In general, an approximate bulk velocity u = (u, v, w) T can be computed from an existing shallow granular model [e.g. [20][21][22] and then a fully coupled model can be developed. An example of how to couple these types of models can be found in Woodhouse et al. [23], where a coupled model to describe the geophysically important granular fingering instability is derived.

Scaling.
Assuming the bulk flow velocity is approximated using such models, the flow quantities are scaled as follows  (4.9) where S r = qL/(HU ) is the non-dimensional number defined as the ratio of the mean segregation velocity, q = g n c s a −ρ γ, to a typical magnitude of the normal bulk velocity.

Solutions for limiting cases.
For the kinematic limiting cases, solutions are constructed for (4.9) for a velocity field u = (u(z), 0, 0) in the domain 0 ≤ z ≤ 1 and x ≥ 0, for a mixture flow of unit height. For simplicity, initially we consider thin/shallow flows in which the velocity profiles are almost linear [24]; hence, to a good approximation we can take the shear rate,γ, to be constant. For the case with S r =S r constant, Eq. 4.9 becomes (4.10) As in Gray and Thornton [18], we use a simplified boundary condition with no flux of particles at the free-surface and the base by considering F (φ) = 0 at z = 0, 1. Solutions to (4.10) can be constructed for two types of inflow boundary conditions, prescribed at x = 0, termed as homogeneous mixture inflow and normally graded mixture inflow. In the former type, a homogeneous mixture of concentration with φ(0, z, t ) = φ 0 , constant, enters at x = 0, whereas in the latter case normally graded particles, i.e., a mixture in which smaller particles lie on top of the larger particles, enters at x = 0.

Analytical solutions.
Assuming the flow has reached its steady state, the differential equation (4.10) is restated as which is a quasi-linear partial differential equation. By using the method of characteristics, as applied in Gray and Thornton [18] and Thornton et al. [19], the above equation (4.11) is analytically solved for specific inflow boundary conditions.

Characteristics.
An exact solution to the above equation is obtained via the method of characteristics, such that φ is constant, φ λ , along each characteristic curve given by The characteristics describe the flux of information into the domain. Given a velocity field, u(z), the above equation is integrable as downslope velocity u is a function of z alone and φ λ is constant. Solutions for general velocity fields can be obtained by introducing a depth-integrated velocity coordinate ψ, where Equation (4.12) can then be integrated to give straight line characteristics, 14) in terms of the mapped variables, with (x λ , ψ λ ) as the starting point for varying φ λ . We choose a scaling such that without loss of generality the mapped coordinate ψ = 1 at the free surface z = 1. As u(z) > 0 is taken from the onset to simplify matters, the map between the physical and depth-integrated coordinates can easily be constructed for a whole class of general velocity fields. Segregation model Jump conditions.
Experimental evidence of segregating flows reveals that concentration jumps or shocks can emerge [15]. The presence of shocks implies that the segregation equation (4.9) or (4.10), is no longer valid because particle concentration φ is then no longer continuous.
Hence, a jump condition should be applied across the discontinuity [25]. We derive the jump condition from an integral version of the conservative form of segregation equation (4.10). We have Assuming a jump in φ exists at z = J (x) and following Whitham [25], the jump condition is with J ′ = d J /d x, '+' and '−' denoting the limits on either side of the discontinuity/jump J (x), and the bracket, '[ ]', denoting the difference of the enclosed function value at upper and lower limit. Since u(z) is continuous, the above equation can be restated as which, when solved, gives the location, z = J (x), of the shock. Following Gray and Thornton [18], we restate the above equation in terms of depth-integrated velocity coordinates, (4.13), i.e., independent of the assumed monotonically increasing velocity profile, with ψ J = ψ(J (x)).

Jumps in mapped coordinates.
For a homogeneous inflow condition and purely size-based segregation, i.e.,ρ = 1, positions of the shocks are determined from the shock relations (4.18). By substituting φ + = 1 and φ − = φ 0 and integrating with boundary condition ψ = 0 at x = 0, the position of the shock is Similarly, by substituting φ + = 0 and φ − = φ 0 into the shock relations and integrating with boundary condition ψ = 1 at x = 0, the position of shock The shock ψ 2 propagates downwards to merge with shock ψ 1 at the triple point, x t r i pl e = (φ 0 + (1 − φ 0 ) s a )(( s a − 1)(g n /c)) and at ψ = φ 0 in depth-integrated velocity variables, resulting in a third shock separating the 100% species-1 and species-2 regions. The shock position, ψ 3 , is determined by substituting φ + = 0 and φ − = 1 into the shock relations which on integrating gives ψ 3 = φ 0 , for x ≥ x t r i pl e . When the shock positions ψ 1 , ψ 2 and ψ 3 are mapped back to physical coordinates, they yield the solid lines in Fig. 4. 1.
Physical solutions and comparison with other models.
The shock positions, in terms of the depth-integrated velocity coordinates, are valid for all velocity fields given that the fields have a single valued map from the depth-integrated In case of a simple shear flow where α = 0, from (4. 19) we find that z = ψ. Thereby, the jump/shock positions are mapped back to the physical space. Fig. 4.1 shows the results from our model, the model of [18] and the model of [15]. To allow direct comparison of the effects of the flux functions, between our model and [18]'s models the segregation velocity,S r = 1, is taken to be unity in both theories. Two things are clear from this comparison: firstly, the point of full segregation (location of the triple point) is further down-stream for our model, by as much as a factor of two for φ 0 = 0.1; secondly, the distance to full segregation is a function of φ 0 in our model (in contrast to the model of [18]). Our flux function F (φ) is convex and more general; it is apparent that the form of this function can have a large effect on the predicted distance to full segregation.

Numerical solutions.
Solutions can be computed to the segregation equation

79
leading to a Bagnold-type velocity profile. Substituting Eq. (4.20), dropping tilde, in the scaled segregation equation (4.9) gives where M = (L/HU )(g n /c)(tan θ − µ c ) 1.5H g cos θ/k(s 1 ) a is a material parameter, µ c is the friction coefficient and k a non-dimensional constant when a = 1 or a dimensional constant when a = 2 or a = 3. We choose M = 0.1 such that the DPM simulation and continuum numerical solution takes the same time to reach steady state. We use a high resolution shock capturing method, the space discontinuous Galerkin finite element method (space-DGFEM) [26], implemented in our open-source DGFEM solver hpGEM [27]. Full details of the package including the numerical codes used to generate Fig. 4 [18] and [19] for flux functions of the form In

No/weak segregation.
One of the key new features of our model is the ability to predict the ratios of size and density for which no or very weak segregation should occur. Due to the non-dimensionalisation used, this is simply given by the line s a = ρ; recall that a = 1 for the pressure function suggested by [13].
To validate our no or zero segregation prediction, we used fully three-dimensional discrete particle simulations (DPMs) (a.k.a Molecular Dynamics (MD) simulations or the Discrete Element Method (DEM) simulations), implemented in our in-house open source DPM package, MercuryDPM. This package has previously been used by Thornton et al. [16] and Weinhart et al. [29] to investigate size-segregation in chute flows. Full details about the MercuryDPM and the source code used for this paper can be found at the website MercuryDPM.org.

DPM simulation Setup
We simulate a homogeneously mixed bi-disperse mixture of particles. We will refer to the two different particles as species type-1 and type-2. In general, if d i , i = 1, 2, . .., is Segregation model defined as the particle diameter of each species type-i , then the mean particle diameter is with φ i being the volume fraction of particles of species type-i . Hence, for a bi-disperse mixture d m = φd 1 + (1 − φ)d 2 with φ being the volume fraction of particles of species type-1. Correspondingly, ρ m is the mean particle density.
In our chosen coordinate system, we consider a cuboidal box, periodic in x and y, inclined at 26 • to the horizontal. The box has the dimensions L × W × H = 20d m × 10d m × 10d m . To create a rough base, we fill the box with a randomly distributed set of particles of both types of diameter d m and simulate until a static layer of about 12 particles thick is produced. Then a slice of particles with centres between z ∈ [9.3, 11]d m are fixed and translated 11 particle diameters downwards to form the base of the box. A solid flat wall is added below this layer of particles to ensure that none of the flow particles will fall through. Once the base is created, the box is inclined and filled with homogeneously mixed bi-disperse mixture, see Fig. 4.5 and [24] for details.
The parameters in our DPM simulations are non-dimensionalised such that mean particle diameter d m = 1, its mass m m = 1 and the magnitude of gravity g = 1. Therefore the non-dimensional mean particle density ρ m = m m / V m = 6/π, with mean particle volume V m = π( d m ) 3 /6 and time scale d m /g , with ' ' indicating non-dimensional quantities. Furthermore, given the diameters and densities of each species type, the particle size and density ratios, s and ρ, is defined as s = d 2 /d 1 and ρ = ρ 2 /ρ 1 , respectively.
The box is filled with a bi-disperse mixture containing particles, where V box = 10×20×10 is the volume of the box. The formulae (4.23) ensures that the ratio of total volume of species-1 over total volume of all the particles is φ and the dimensionless height of the flow, H is the same for all simulations. Using (4.23), for homogeneous initial conditions (randomly mixed), with initial particle volume fraction φ i = 0.5, a series of DPM simulations for different values ofρ andŝ were carried out to Segregation model For all the performed simulations, we use a linear spring dashpot model with a contact duration of t c = 0.005 d m /g , coefficient of restitution r c = 0.88 and contact friction coefficient µ c = 0. 5

Analysis
The sensitivity to both basal and initial conditions on the steady-state has been thoroughly investigated and hardly any sensitivity was found [32]. Once the flow has reached its steady state, we calculated a relative difference between the centres of mass D com (ρ,ŝ) = (COM2 − COMB)/COMB, as a function ofŝ andρ. Here COM2 is the vertical centre of mass of species-2 particles and COMB is the bulks' vertical centre of mass. For a given s andρ, the flow is steady when the function value, D com , remains constant with time.
In Fig. 4.4, we plot the values of D com for givenρ andŝ. When the value of D com is positive, particles of species-2 are near the free surface; vice-versa, when it is negative, particles of species-2 are near the base, see Fig. 4. 3. Close inspection of the data shows very weak segregation along the solid line s a = ρ with a = 3, also implying that the pressure is scaled by the volume of the particle. Below the solid line the species-2 particles rise towards the free-surface and above the solid-line species-2 particles fall towards the base. The dashed-solid line corresponds to the prediction via kinetic theory for a binary mixture [33]. The mismatch between the dashed-solid line and the solid line is probably because the prediction by Jenkins and Yoon [33] is valid only for mixtures in which large particles are dilute in a dense mixture comprised of smaller particles. This assumption is best satisfied for lowŝ andρ limit; which is clearly seen in the bottom left corner of Fig. 4. 4. In this limit, our continuum theory is not valid, as from Savage and Lun [15] and Thornton et al. [16] it is known that for size ratios greater than 2 percolation effects are present. This corresponds to s < 0.5, in our current study, and s > 2. Percolation of small sized particles through the matrix of large sized particles occur simply as a result of gravity, in absence of shear. In contrast to kinetic sieving that requires the combination of gravity and shear to be acitve.
Felix and Thomas [11] experimentally investigated density and size segregation for both chute flows and rotating drums. It is not possible to directly compare their chute flow results with our DPM simulations, as in their chute experiments the type-2 species' (tracers) volume fraction is 10%. It is very hard to directly fit their experimental data to obtain the value of a; however, several qualitative similarities can be found between our DPM simulations and their experimental findings. Firstly, for ( s = 1, ρ ≈ 1.1) they conducted two experiments for which they observed segregation only in one and homogeneous mixture state in the other, indicating that the mixture weakly segregates at that ratio of size and density. For these ratios, our DPM data shows weak segregation (see Fig. 4.4), i.e., the flowing mixture only partially segregates with a fairly homogeneous mixture still existing in the bulk (see Fig. 4.5a). We only see almost complete segregation in the regions of Fig. 4.4, with dark blue or red spots; however, as we measure the centre of mass of both species we can detect even weak segregation, compare the left and right panels of Fig. 4.5 that correspond to cases where D com = −0.1 and 0.5, respectively. The weak segregation in our phase plot seems to correspond to the homogeneous states of Felix and Thomas [11], as such weak segregation could not be captured by the experimental techniques they employed. Secondly, for the cases where ρ ≈ 1.3 and s = (1.0, 1.1, 1.2, 1.3), they observed undeterminable/homogeneous mixture states, which from our DPM simulations again correspond to a weakly, not fully, segregated region in our phase plot. However, for larger size ratios, i.e. s > 2, an opposite trend versus density ratio is observed in their experiments indicating that the segregation reverses direction. Given these differences, we suggest that this behaviour could be due to percolation effects and different volume fractions considered; however, further study is required to confirm this.
Moreover, a generalised pressure scaling function could be obtained directly from the DPMs. This work is ongoing and some early results can be found in [29]. In this paper, we have not incorporated the diffusive nature of these segregating flows into our model [34] and therefore the model predicts full segregation, eventually, on either side of the no segregation line. However, if we had included diffusion in the model then D com could be reinterpreted as a segregation Péclet number (ratio of segregation to diffusive strength) as in Thornton et al. [16].

Summary and conclusions.
In this paper, we focus on bi-disperse flows over inclined channels with size ratios less than 1. 5. We derived a generic continuum model to predict the extent of segregation in a bidisperse granular mixture flow due to differences in particle size and density. For a given density and size ratio, the model predicts the extent of de-mixing in gravity-driven chute flows. For purely size-based segregation, the model has been compared to two previous models. The derived model, for constant shear, is solved analytically using the method of characteristics, and numerically using a discontinuous Galerkin finite element method. The model was also used to predict the ratios of particle size and density for which very weak/no segregation would occur. This prediction is independent of the details of the drag coefficient between the particles and the bulk velocity profile of the flow. To validate this prediction, we performed discrete particle simulations as an alternative to laboratory experiments of field measurements. The model performs surprisingly well, when compared with the discrete particle simulations, with the fitting parameter 'a' determined from the DPMs. The advantage of our continuum model is that it permits analytical and fast numerical solutions. However, for size ratios s < 0.5 and s > 2, the current model is not valid as percolation effects [15,16] are not considered.

In Chapter 4, we have seen a beautiful blend between a mixture theory continuum model and discrete particle simulations, where the two theories were used to model and understand bidisperse granular mixtures flowing over rough inclined channels. However, despite the illustration of an efficient micro-macro transition technique for multi-component mixtures (chapter 3), a direct application of this mapping tool is still something we have not yet considered in the previous chapter(s). Thereby, this chapter focusses on utilising the micro-macro transition technique.
To be submitted.

Introduction
Granular mixtures often tend to arrange themselves in certain patterns when stirred, shaken or sheared [1]. These mixtures often comprise of constituents varying in size, density, inelasticity, shape, surface roughness, etc. When such polydisperse mixtures are subjected to external forces, individual studies [2][3][4][5][6] confirm the influence of these constituent properties in forming patterns, e.g., segregation or de-mixing. However, in free-surface flows over inclined channels differences in size and density are the primary factors for segregation. Among several mechanisms responsible for segregation, kinetic sieving [7] is the dominant one in dense granular flows. Kinetic sieving is a trivial mechanism. As the mixture flows down the inclined channel, fluctuations in the local pore space cause smaller particles to fall into the space/gaps created beneath them. The fine ones, i.e. small sized particles, easily fit into these pores leading to gradual percolation of them towards the base of the flow. Simultaneously, force imbalances lever/squeeze the larger particles towards the surface. This simple mechanism results in stratified layers that one terms as segregation, as seen in Chapter 4. Opposing kinetic sieving is diffusive remixing, which causes random motions of particles as they collide and shear over each other [8]. Based on the relative strength of either mechanisms, the mixture strongly or weakly segregates. Apart from kinetic sieving, which is a purely size-based effect, buoyancy effects due to differences in particle density also play a major role in particle segregation [9]. For bidisperse mixtures, varying in particle size and density, experiments [10] and numerical simulations [11] (Chapter 4) indicate a balance between the two driving mechanisms, i.e kinetic sieving and buoyancy effects, which in turn keeps the mixture homogeneously mixed. However, a detailed understanding of the dynamics of segregation becomes more difficult with the increase in the number of particle properties, i.e. size, density, shape etc., one includes in a study.
To get a grip on the segregation dynamics, as an alternative to experiments, we employ discrete particle simulations (DPMs) [12]. From DPMs, macroscopic quantities corresponding to the macroscopic field variables, appearing in continuum-mechanical models, are extracted. To do this, we use an appropriate micro-macro transition procedure called coarse graining, which has been described in detail in Chapter 3 and references therein ( e.g. [13,14]). This mapping technique has been successfully applied to flows near boundaries or discontinuities [15], shallow granular flows [16] and bidisperse mixtures with varying particle size [17]. By this micro-macro mapping, one can easily use the DPMs as a calibration tool for the continuum models.
Based on the existing knowledge pool associated with quantifying segregation dynamics, a few continuum models [11,[18][19][20][21]23] have been formulated to predict segregation in gravity-driven granular flows. Out of these, we focus on mixture theory [24] based kinetic sieving models [21,25]. Recently, the Gray and Thornton [21] model was subjected to verification by extending the coarse graining technique [15] to mixture flows [17]. As a result of this, for a certain particle size-ratio, Weinhart et al. [17] observed that (i) the large particles do support a fraction of the stress larger than their volume fraction; (ii) a simple linear interspecies' drag law is not sufficient to describe bidisperse mixture flow dynamics; (iii) the smaller particles support a fraction of the kinetic stress larger than their volume fraction, as postulated by Fan and Hill [22] in their shear-induced seg- Figure 5.1: A snapshot of a bidisperse mixture flowing in a periodic box inclined at 26 • to the horizontal (discrete particle simulation). Colours/shades indicate the base/boundary (yellowish green, F b ), species type-1 and type-2 (blue, F 1 and red, F 2 ). We define the bulk as F 1 ∪ F 2 . regation theory.

Simulation setup
In this chapter, we extend the purely size-based Gray and Thornton [21] model to consider density differences as well, the derivation steps are similar to the steps followed in Chapter 4, however, we employ a simpler pressure scaling function. Similar to what was done in Weinhart et al. [17], on using the coarse graining expressions from Chapter 3, we verify if the observations of Weinhart et al. [17] are still valid, when considering bidisperse mixtures varying both in size and density. Furthermore, from the discrete particle simulations, we also determine the closure parameter, which is dependent on both the size-and density-ratio.

Simulation setup
To illustrate the micro-macro technique, as an example, fully three-dimensional (3D) discrete particle simulations (DPM) are used as an alternative to experiments, to investigate segregation dynamics in bi-dispersed mixture flows over inclined channels. For setting up the simulations, we use our in-house open source DPM package, MercuryDPM. Previously, this simulation package has been used by Thornton et al. [26] and Weinhart et al. [17] to investigate size-based segregation in inclined channel flows.
We consider a cuboidal box, periodic in x-and y-direction, inclined at 26 • to the horizontal. The box has dimensions L × W × H = 20d m × 10d m × 10d m . For a box to have a rough base (bottom), we fill the box with a randomly distributed set of particles with diameter d m and simulate them until a static layer of about 12 particles thick is produced. Then a slice of particles with centres between z ∈ [9.3, 11]d m are fixed and translated 11 particle diameters downwards to form the rough base of the box. To ensure no flowing particles fall through the base, a solid wall is placed underneath this static layer. Once the rough base is created, the box is inclined and filled with a homogeneously mixed bi-disperse mixture of particle diameters d 1 and d 2 , and densities ρ 1 and ρ 2 , as illustrated in Fig. 5.1, see [16] for more details.
In our DPM simulations, we non-dimensionalise the parameters such that the mean particle diameter d m = 1, the mean particle mass m m = 1, the magnitude of gravity g = 1. Implying, that the mean particle density ρ m = m m / V m = 6/π and the mean particle volume V m = π( d m ) 3 /6. The non-dimensional quantities are denoted by ' ' . Furthermore, given the diameters and densities of each particle species type, the particle size and density ratio is defined as s = d 2 /d 1 and ρ = ρ 2 /ρ 1 , respectively.
Finally, we fill the box with the bi-disperse mixture comprising

Mixture theory
We choose a coordinate system such that the x-and y-axis points in the downslope and the cross-slope direction of the channel, and the z-axis points in the upward direction normal to the channel, see Fig. 5. 1. Similar to Chapter 4, we briefly introduce the mixture theory postulates. For a bidisperse mixture of grains, mixture theory [24] for a binary continuum postulates that both constituents of the mixture simultaneously occupy both space and time. This results in overlapping fields leading to partial flow quantities, such as partial stress σ ν , density, ρ ν , and velocity, u ν = [u ν , v ν , w ν ] T in the three coordinate directions, corresponding to each constituent indexed ν = 1, 2. Each of the constituents satisfies the following fundamental balance laws of mass and momentum stated in terms of the partial fields as where g = (g t , 0, −g n ) T is the gravity vector, with g being the standard acceleration due to free fall; g t = g sin θ and g n = g cos θ. The variable β ν represents the interspecies drag force due to resisting motion between the constituents. As these forces are internal, from Newtons' third law the sum of these drags must be zero, i.e., β 1 + β 2 = 0.
Given a unit mixture volume, including the interstitial pore-space, each of the constituents occupies a volume fraction φ 1 or φ 2 . Hence, by definition, the individual volume fractions sum to unity φ 1 + φ 2 = 1. Furthermore, the bulk density ρ, barycentric granular velocity or bulk velocity u and the bulk stress σ are defined as ρ = ρ 1 + ρ 2 , u = (ρ 1 u 1 +ρ 2 u 2 )/ρ , σ = σ 1 +σ 2 respectively. The intrinsic variables defined for each of the constituents do play an integral role in the constitutive theory. These quantities are related to the partial quantities and hence are the key features of the mixture theory. The intrinsic density of the constituent ρ ν * , i.e. the mass of the constituent per unit volume of the constituent, is related to the partial density by constituent volume fraction φ ν . The same applies to the partial pressures of the constituents. However, in standard mixture theory the partial velocity of the constituent is identical to the intrinsic velocity of the constituent where all the intrinsic quantities are denoted with ' * '.

Gravity-driven segregation
Most of the mixture flows involved in industry and geological applications are shallow in nature, implying that the flow quantities in the downslope and cross-slope direction are nearly uniform. Moreover, we assume that the partial densities and momenta become quasi-steady even before the flow segregates, implying that the temporal derivatives ∂ t (ρ ν ) and ∂ t (ρ ν u ν ) vanish after a certain equilibrium time t e . On employing the shallowness argument, as done in the previous chapters, from the momentum equation (5.2) 2 , we arrive at Summing (5.4) over each particle species-type ν = 1, 2 for α = z, setting σ zz | z=inf implies that the flow is in lithostatic balance ∂σ zz ∂z = −ρg cos θ. (5.5) The key idea behind the kinetic sieving model of Gray and Thornton [21] is that as the small particles percolate downwards through the granular assembly, they carry less of the weight and the large particles carry proportionately more. To quantify this, they introduced a new scaling function f ν , which determines the amount of overburden pressure (weight) is to be distributed among each of the constituents. In standard mixture theory, the constituent pressure is assumed to be linearly related to the bulk pressure through the volume fraction, i.e. f ν is assumed to be equal to φ ν . Here, we have, however, a crucial deviation from the standard approach in order to account for the effects of segregation, i.e. we define the intrinsic partial pressure as The stress fraction f ν is defined using the following functional form with φ 1 f 1 + φ 2 f 2 = 1, f 1 = 1 when φ 1 = 1 and f 1 = 0 when φ 2 = 1. Similarly, f 2 = 1 when φ 2 = 1 and f 2 = 0 when φ 1 = 1 In order to use these scalings, one still needs to determine the magnitude and sign of B ν . For a bidisperse mixture varying in size alone, Gray and Thornton [21] considered the pressure scalings to be f s = 1 − Bφ l and f l = 1 + Bφ s , (5.8) where superscripts s, l denote small and large constituents, with B s = −B and B l = B. Furthermore, Gray and Thornton [21] also considered the interaction drag or interspecies' friction to take the form of Darcy's law, where c is an inter-constituent drag coefficient. The inter-particle surface interaction force is given by p∇(φ ν f ν ), thus ensuring that the particle percolation is driven by the partial pressure gradients. Substituting the expressions for the drag force (5.9) and lithostatic balance (5.5) into the normal momentum balance equation for constituent type-ν (5.4) -neglecting normal acceleration terms -results in the species percolation velocities For ν = 1, the above expression is further simplified, for details see Appendix D.2, to By substituting the above percolation velocity (5.11), with u ν = u and v ν = v (see Appendix C), in the continuity equation for ν = 1 leads us to the segregation governing equation 12) where ρ = ρ * 2 /ρ * 1 is the particle density ratio. For a bidisperse mixture varying in size alone, ρ = 1, we arrive at the purely size-based segregation model of Gray and Thornton [21]. However, in order to utilise the above segregation model one still needs to determine the values of B and c. In this chapter, we determine the values of B using discrete particle simulations (Sec. 5.2) and the coarse graining expressions from Chapter 3,

Results
At around t e ≈ 3.1s, all the flows (DPM simulations) considered in this chapter become quasi-steady implying that the flow density, velocity and pressure profiles change at a very slow rate. Once the mixtures flow over a duration of t e , particle species of type-1 and type-2 either rise towards the surface or fall towards the base of the flow, depending on the particle size and density ratio. After around t ≈ 20s, the process of demixing or segregation is complete and steady. With the flows in steady state, the corresponding data sets are coarse grained to obtain the macroscopic flow fields.
In order to carry out the required micro-macro transition, i.e. to determine the value of B, we require the coarse graining expressions for the partial stresses (pressure) and  Figure 5.2: For a bidisperse mixture with s = 1.4 and ρ = 1, we plot the stress fractions φ ν f ν con and φ ν f ν ki n as a function of φ ν and fit it to (5.8). The above fits correspond to B 1 con = −0.005 and B 1 ki n = 0. 40. volume fraction, as the ratio of the partial to bulk pressure is where σ ν zz is total partial stress in the normal depth-direction. We use the word total because, in Chapter 3, we defined the total partial stress (3.39) as the sum of both partial contact (3.37) and kinetic stress (3.40) fields, σ ν zz = σ con,ν zz + σ kin,ν zz . Using the coarse graining expressions illustrated in Chapter 5, depth-profiles of volume fraction and partial stresses are constructed by averaging both spatially (x-and y-direction) and temporally. Utilising these fields and substituting (5.8) in (5.13), the sign and magnitude of B in the scaling function f ν is simply given by Furthermore, in order to see whether segregation is primarily caused by the contact or the kinetic part of the total stress, we compute the contact and kinetic stress fraction To illustrate this, we consider three bidisperse mixtures with different particle size-and density-ratios, ( s, ρ) = { (1.4, 1.0), (1.0, 1.4), (0.5, 0.6) }. In Figs. 5.2 and 5.3, we plot the stress fractions φ ν f ν con and φ ν f ν kin against the volume fraction φ ν for bidisperse mixture flows varying both in size or density. It follows, that for bidisperse mixtures varying in size alone, the smaller particles (type-1) have a kinetic stress fraction slightly higher than their volume fraction. This is consistent with results in [17] and in agreement with shear induced segregation theory of Fan and Hill [22]. Similarly, for bidisperse mixtures varying in density alone, the heavier particles have a kinetic stress fraction higher than their volume fraction. For a bidisperse mixture varying in both size and density, the behaviour of the stress fractions against volume fraction is illustrated in Fig. 5. 4. It follows, that the smaller and heavier particles (type-1) have a kinetic stress fraction slightly larger than the volume fraction.

Summary
We illustrated a simple application of a micro-macro transition using the coarse graining expressions of Chapter 3. By employing these we not only determined the closure parameter, B ν , for our formulated continuum model, but also showed that by utilising this mapping technique one can better understand the segregation dynamics or mechanisms.
For bidisperse mixtures varying in both size and density (i) the smaller and heavier particles do support a fraction of the stress (kinetic) larger than their volume fraction; (ii) a simple linear interspecies drag law is not sufficient to describe bidisperse mixture flow dynamics. There exist, however, many more aspects that still need to be better understood regarding particle segregation, such as determining the correct pressure scalings, diffusive remixing and many more, which will be the subject of our future studies. 97 Figure 6.1: Overview of the monotonically increasing research interests concerning discrete particle simulations of particulate systems in the last 10 years. The statistics are obtained from Google scholar using the following keywords: 'discrete element method', or 'discrete particle simulation'. Of course there are other keywords that could be used such as 'discrete particle method', however, the main intention was to illustrate the rapid increase in utilising particle simulations.

Background and Aims
Granular materials -assemblies of multiple, discrete particles or 'grains' -are, after water, the second most manipulated material on the planet [1], playing important rôles in multitudinous natural and industrial processes [2][3][4][5][6][7][8][9]. Yet in spite of their ubiquity and industrial relevance, and despite the significant volume of research in the field [10], the behaviours of granular materials remain poorly understood and difficult to predict [11]. Our inability to accurately predict and control the behaviour of granular flows carries many negative consequences for industries required to handle particulate media. For instance, phenomena such as jamming, agglomerating or clogging [12] may significantly reduce the efficiency of industrial flows [13] or even cause catastrophic structural failures in storage silos [14,15]; phenomena such as segregation (gradual separation of dissimilar particle species) [16,17], meanwhile, can prove similarly problematic, for instance in the pharmaceutical industry [18], where a lack of mixing may -for obvious reasons -have serious repercussions.
In order to avoid issues such as those described above, companies often resort to ad hoc solutions, such as the construction of scaled down test systems or pilot plants, or even opting to 'hit and hope' -simply constructing a system and attempting to resolve any issues later. The former solution can often prove extremely costly and, since the scaling behaviour of granular materials is also incompletely understood [3,19,20], may still not provide any information relevant to the full-scale set-up. The latter, meanwhile, is likely to result in decidedly sub-optimal efficiency or, worse, might fail entirely.
Increasingly, industries are turning away from such 'practical' methods, choosing instead to utilise computational models. In particular, simulations produced via the discrete particle method (DPM) [21][22][23][24][25][26][27][28][29] -or discrete element method (DEM) as it is more commonly referred to -are used in order to gain an understanding and predictive, quan-titative knowledge of the dynamical and physical properties of granular systems. Discrete particle models -a rapidly evolving field of research, see Fig.6.1 -carry a significant advantage over practical investigations (experiments) in that they can provide a plethora of information regarding the dynamics of each particle within the system at any given point in time. This level of detail is, at present, near-impossible to achieve in practical, experimental settings -techniques that can capture information regarding all particles within a system typically suffer from poor time resolution [30], while systems capable of imaging in real-time can generally only image relatively small systems, or limited fractions of larger systems [31][32][33]. Similarly, techniques which provide detailed information regarding the dynamics of a system [34] are often limited in their ability to provide information regarding detailed structural and contact properties, and vice-versa.
Unfortunately, discrete element methods also carry certain notable disadvantages. Firstly, the simulation of systems containing large ( 10 6 ) numbers of particles will typically require hardware that may be prohibitively expensive and timescales that are often unfeasible. However, despite the current trends of increasing processing power in computational systems and the falling cost of technology, such a problem is likely to persist for a considerable period of time. This being the case, a more time-efficient alternative is to use discrete element methods only when necessary, i.e. to determine the closure or constitutive relations necessary for a continuum approach (micro-macro mapping -see [35,36] and references therein). Alternatively, it is possible to couple simulations with continuum methods (coupled multi-scale methods, CMSM -see [37][38][39]).
An additional disadvantage of the discrete element method is the necessity to implement certain modelling assumptions when emulating particulate systems, e.g. particle shapes (spherical or non-spherical), type of interactions (contact models) etc., which are often approximations compared to 'real life'. This raises an important question:

How accurately can DEM/DPM models emulate 'real' granular systems?
It is this question that forms the central theme of this paper. The ability of DEM simulations to faithfully recreate the behaviour of physical systems is pivotal to their use as viable predictive models. While we have focussed thus far on the relevance of DEM simulations in industrial settings in order to quickly and simply give a feel for the magnitude of the issues at hand, computational models also have numerous other highly important applications regarding questions of science and progress in other disciplines. For instance, granular systems provide useful, easily accessible analogues for various biological systems [40,41], which may be difficult to study in-vivo or in the field, allowing us to extract valuable information regarding, for instance, the behaviour of bacteria [42], animals [43,44] and even humans [45][46][47]. Computer models of granular systems may even eventually provide us with the ability to predict and prevent disastrous geophysical events, such as earthquakes, avalanches and landslides [48][49][50].
Before any of the above can be attempted, we must first ensure that computer (particle) simulations can indeed provide an accurate representation of the real systems that they are designed to represent. In this work we raise, and aim to resolve, a series of questions pertaining to this matter: what is the current state of the field of discrete element modelling? Under what conditions, and for what parameter ranges, can specific sys-tems be accurately modelled at present? How do we ensure that we implement sensible values for particle properties? How can we test the validity of the numerous simplifying assumptions used in simulations? Through an analysis of the existing literature and with reference to our own experimental and numerical work we aim to address these vital and pressing questions.

Article Outline
This chapter is structured as follows: In the following (Sec. 6.2), we begin by briefly listing recent developments relating to the methods used to simulate particulate systems, with particular emphasis on the post processing techniques implemented. We also detail current works pertaining to parallelised discrete element simulations (Sec. 6.2.1) and micro-macro mapping techniques (Sec. 6.2.2). We then discuss how discrete element simulations are used in conjunction with experimental systems and the data produced by these systems (Sec. 6.3), detailing the typical process by which suitable quantities can be determined to recreate the real-world system being modelled. We then assess how well DEM results compare to experimental data for two common experimental setupsvibrated beds (Sec. 6.3.1) and inclined channels (Sec. 6.3.2) -and for a variety of system parameters. Finally, we summarise our results and provide an outlook (Sec. 6.4).

Particle simulations
Over the past decades, a range of computational methods have been developed to address and understand granular dynamics [51]. These include cellular automata (CA), direct simulation Monte Carlo (DSMC), contact dynamics (CD) and the discrete element/particle method (DEM/DPM).
In a simple, deterministic technique such as lattice-based cellular automata (CA), particle positions are determined using rule-based mathematics or physics-based equations. CA utilise lattice like structures (grids) in which the physical domain is divided into cells. Each cell corresponds to one of the defined number of states (see Goles et al. [52]). The technique has been used, previously, to understand the flow of granular materials in silos [53][54][55], sand piles [51,52,56], and annular shear cells [57,58], as well as the granular flow over inclined channels [59,60], in rotating drums [61,62], and in hopper flows with irregular particle shapes [63]. Although the method is computationally fast, the predicted particle dynamics lack quantitative accuracy, agreeing only qualitatively with experimental observations. Granular dynamics have also been modelled using a stochastic direct simulation Monte Carlo (DSMC) approach. This method solves the Boltzmann equation based on a direct statistical simulation of the molecular processes described by the kinetic theory, see [64,65]. The technique, first proposed by G. A. Bird [64], is both computationally efficient and relatively simple to implement as compared to DEM implementations. The method has been successfully utilised to simulate dilute granular flows in two-and three-dimensional vibrated containers [66] and Couette flows [67], for instance. Recently, modified DSMC algorithms have been developed [68][69][70] which, due to their improved algorithms, are faster than traditional DSMC methods and, more importantly, consistent with DEM results [70].

Discrete element method (DEM)
With the advent of computing technology, research in granular media has increased fourfold since the past few decades. It was in the late 1970s', when the myth of tracing the trajectory of a particle or grain metamorphosed into a much desired reality. Since then, particle tracking has become essential in numerous applications involving particulate media.
DEM is categorised into two main -hard-and soft-particle [21,51,71] -categories. In the hard-particle (rigid particles) approach, collisions are assumed to be instantaneous and binary whereas in the soft-particle approach the collisions are considered to be enduring, i.e. the contact duration is finite 1 . Hard-particle based computations are, in many cases, significantly faster as compared to equivalent soft-particle simulations. Even today, when necessary, studies still utilise the hard-particle based particle simulations, e.g. in dilute, rapid collisional flows. However, for highly dense, quasi-static flows, where particle contacts are enduring and multiple, the soft-particle approach is, clearly, preferable. For a detailed review of these approaches, their corresponding algorithms, theoretical advances and applications we refer to recent reviews [21,51,[72][73][74] and the references therein.
Even though the majority of the references listed in these articles [21,51,[72][73][74] consider spherically shaped particles, particulate systems encountered in industry or nature are often made up of non-spherical and/or irregularly shaped particles. Studies have shown that particle geometry considerably affects the bulk dynamics of particulate mixtures. As a result, a variety of techniques have been formulated to construct different particle shapes (ellipses/ellipsoids, super-quadrics, polygons/polyhedrons, multispheres etc.) and to perform efficient contact detection for these more complex geometries. However, we do not review these state of the art techniques as these have already been comprehensively addressed [71,72,[75][76][77][78][79][80].

Computational speed-up
Despite tremendous efforts to allow DEM to successfully emulate reality, the computational time associated with soft-sphere DEM has always been considered a liability. While soft-sphere based DEM provides us with a plethora of useful particle data, this data requires a considerable amount of time to compute. For example, let us consider a DEM simulation of a granular mixture. During a dynamic simulation process, one has to store and update the particle positions at every time step, and also needs to detect all collisions or contacts between particles and calculate the resulting forces due to these contacts. However, in case of static mixtures, e.g. simulations concerning stationary assemblies of rocks, it is the force computations that consume more time since the detection of contacts for fixed neighbours is just a one time operation. Nevertheless, both contact detection and force computations carry a high computational cost in most cases.
As a stepping stone towards minimising this computational effort, recent studies have succeeded in developing efficient solutions (algorithms) either by developing state of the art contact detection algorithms [81][82][83][84] or utilising multi-core [85][86][87] processors, heterogeneous CPU 2 -GPU 3 architectures [87][88][89][90], supercomputers [91] or parallelised GPU clusters [92]. A modern day personal desktop -possessing two or four cores -is able to simulate a fully three-dimensional simple soft-sphere (linear force model) DEM simulation comprising 10,000 to 200,000 particles within a couple of hours. A highly parallelised cluster [93][94][95][96], meanwhile, would only require a couple of minutes to perform such a task. Despite of advances in parallelised DEM packages, simulation of current industrial applications which involve multitudes of highly polydisperse particles, in the order of 10 9 , is difficult because of the inherited complexity in the detection of contacts [86]. As a result, studies have been focussing on a GPU-based framework for developing a highly parallelised GPU based DEM solver. See [92] for a brief performance overview for CPU-and GPU-based systems.
With such massive parallelism now available using programmable GPU hardware, rapid advancements have been witnessed during the past 10 years. Innovatory simulations on the order of millions of particles are being simulated (conducted) in quasi real-time in confined environments [97], rotating drums [98], blenders [99] and granular soils [100], to give but a few examples. This is thanks to the sophisticated algorithms developed for efficient collision detection [84,[101][102][103][104][105][106][107] and the construction of memory-efficient data structure [108]. Using high performance GPU, [109] investigated the effects of differing particle sizes in granular mixture flows, whereas [110] used it to simulate fractures in heterogeneous media. Additionally, to bridge the gap between ideal and realistic mixtures, studies have also considered the simulation of non-spherical particles using the GPU-based framework, see [111] for triangular particles, or [112,113] for convex polyhedrals.

Micro-Macro transition
Given that DEM is an efficient tool utilised to probe the intricate details of granular dynamics on the micro (particle and contact) scale, macroscopic quantities, such as density, velocity, stresses and other necessary fields, are essential in any study involving development, validation and/or calibration of a continuum model. Besides using this mapping tool to calibrate or validate continuum models, macroscopic fields are useful to quantitatively compare experiments and particle simulations as well. For example, let us consider monodisperse flows (mixtures made up of particles of the same type) over inclined channels. In such experiments, often techniques such as particle image velocimetry [115] are used to obtain velocity fields [116]. In order to compare these fields with the DEM simulations, a micro-macro technique is essential.
Mapping of the microscopic scale information onto a macroscopic scale, e.g. continuum, has been under focus since the classical studies by [117,118] and others [119]. Based on a variety of theoretical postulates, various methods for local averaging have been formulated to extract these macroscopic quantities efficiently, for instance the binning of the microscopic fields into small volumes [73] and the method of planes [120]. However, most methods are restricted in terms of their application due to various limitations, (see [73] and the references therein). One of the requirements for multi-scale methods is to efficiently map the particle dynamics (microscopic) onto a macroscopic 2 Central processing unit 3 Graphics processing unit Figure 6.2: Comprehensive illustration of the contemporary philosophy adopted by many studies concerning utilisation, calibration and validation of the discrete element method [114]. (a) Micro-macro transition (b) Two-way coupling of the continuum and discrete model, also defined as coupled multi-scale mechanics (CMSM). field, which in turn satisfy the classical equations of continuum mechanics, i.e. the fundamental balance law of mass and momentum D t (ρ) + ρ∇ · u = 0, D t (ρu) + ρu∇ · u = ∇ · σ + ρg, (6.1) where D t is the material derivative. The above equations are stated in terms of mass density ρ, bulk velocity u, and the stress tensor σ. "Coarse graining" 4 or "homogenisation" approaches to granular materials first appeared in the work of [121] and have been extended in a number of studies [122][123][124][125][126][127][128][129][130][131][132][133]. Coarse graining techniques have two essential advantages over other types of averaging techniques. Firstly, the macroscopic quantities produced exactly satisfy the continuum laws of motion. Secondly, they are applicable to both static and dynamic granular media. With these advantages, the coarse graining approach has been utilised to study the results of computations or experiments and characterise them in terms of their density, velocity, stress, strain, couple-stress and other fields [134][135][136][137][138][139][140]. Futhermore, the coarse graining method described in [133] has been extended to granular mixture flows near boundaries or discontinuities [141,142] and bidisperse mixtures [36,143]. These extensions [141,143] have been applied to analyse layered flows [35,144] and segregation phenomena in bidisperse granular mixtures [36].

Discrete Particle Simulations and Experimental Systems
Even with the innovative theoretical extensions discussed in Sections 6.2.1 and 6.2.2, emulations of 'real' particulate systems or experiments is still non-trivial. Fig. 6.2 illustrates two conceivable modelling philosophies that utilise discrete element simulations. In Fig. 6.2a, we see a conventional approach adopted by the majority of contemporary studies dealing with validation and calibration of discrete element simulations using experiments. Single particle experiments are vital for establishing the contact parameters (material properties) essential in discrete element simulations. Once appropriate values for the contact parameters are established, both small and large scale simulations are set up. Through small scale simulations, one is able to understand the intricate dynamics -by micro-macro mapping -essential for developing an accurate predictive continuum model. Once the continuum model is formulated, it can then be validated using a large scale simulation or experiments. There exist several granular continuum formulations that need closure or constitutive relations which are unknown, e.g. [145][146][147], and need to be determined experimentally or via DEM simulations, e.g. [35,36,148]. However, when modelling more complex particulate systems, a two-way coupled multiscale mechanics (CMSM) approach is adopted, see Fig. 6.2b. Single particle experiments are utilised to determine the particles' material properties, which in turn are incorporated into small scale discrete element simulations. With the right material parameters, the two simulation models (DEM and continuum) are dynamically coupled such that a two-way feedback exists between the two models. The idea is to use a macroscopic continuum model even where the continuum method fails. In this case, the microscopic model, i.e. DEM, solves for the non-continuum part locally and constructs meaningful macroscopic data, which is incorporated into the continuum model. The coupling is performed in select regions in space and time, thus reducing the computational expense and allowing for the simulation of complex particulate systems. A detailed review of this method is given by Weinan et al. [37].
Continuum theories are advantageous over purely simulational methods for several reasons: (i) they describe the processes occurring in particulate mixtures in terms of differentiable quantities, thus allowing the computation of the solution to the problem by methods of mathematical analysis, e.g. finite difference schemes, finite element methods, finite volume, etc.; (ii) they overcome the need to specify exact microscopic configurations -the microscopic effects are captured in terms of macroscopic coefficients or closure relations, e.g. [35,149]; (iii) they can be scale independent of the particle number; (iv) they do not require an in-depth knowledge of DPMs or experiments. Therefore, in order to have accurate continuum formulations, which are more time-efficient, it is essential that the discrete element simulations emulate -qualitatively and quantitatively -realistic phenomena observed experimentally -thus, bringing us to the following section.

Comparing Experiment and Simulation
If discrete element simulation results are to be of any value to research or industry, they must first be validated through comparison with data pertaining to corresponding 'realworld' systems. Such validation can be achieved, to varying degrees of reliability, by assessing the degree to which they compare to one or more of the multitudinous quantities that characterise the state and dynamics of a granular system.
In the following sub-sections, we present examples from various studies in which the behaviour of simulated and 'real' granulates are directly compared, focussing on two widely-researched experimental set-ups, which are commonly employed in industrial settings -vibrated beds and chute flows. By making such comparisons, we aim to build up a picture of those situations in which current simulation models can be 'trusted' to faithfully reproduce real-world systems and, conversely, cases in which results produced by such models are likely to be unreliable.

Vibrated Systems
Vibrated systems are commonly used in industry for various processes, such as the compaction [150], sorting and mixing [151] of powders and granulates, but can also be exploited to provide useful 'model systems' for various natural and biological processes [152,153].
The behaviour of these systems is largely governed by the details of the vibrations used to energise the system -frequency, f , amplitude, A and hence the dimensionless , and the dimensionless energy input parameter S = where N l is the number of layers in the system and ε the normal coefficient of restitution.
These systems are highly appropriate testing grounds for validating DEM simulations for two main reasons: firstly, compared to 'real-world' systems, and indeed other simplified flow geometries where complex surface and frictional effects are more dominant relative to collisional interactions, vibrated beds are relatively simple. Secondly, experimental set-ups containing relatively few particles can be made to exhibit most of the important phenomena observed in larger-scale systems [154,155], allowing simulations to be performed quickly and easily. Vibrated systems can also -depending on the relevant energy input and dissipation control parameters [156,157] -achieve numerous, vastly differing physical states, ranging from jammed, solid-like and glassy configurations exhibiting slow dynamics [158,159] to dilute, highly energetic gases [160,161], thus allowing agreement between experiment and simulation to be verified under various conditions.
We begin by discussing the relatively simple case of dilute, monodisperse systems of spherical particles, where the majority of collisional interactions are binary and contacts between colliding particles are not enduring [162]. In such systems, the specific details of a particle's frictional and elastic properties are comparatively less important [163,164] and, hence, one may expect simulations to be able to recreate the behaviour of experimental systems with comparative ease. Work by Géminard and Laroche [165] looks at what may be considered the 'extremal case' of a dilute granulate -a single particle bouncing on a vibrating plate. Their research demonstrates, for this fundamental situation, a strong, quantitative agreement between experimental and simulational results concerning the mean energy possessed by a particle and its scaling with the relevant driving parameters. Specifically, agreement was observed so long as the dimensionless acceleration, Γ, with which the system is driven is greater than or approximately equal to 1.4 (the situation for which, in this system, particle motion is no longer synchronised with base motion). Here, the dimensionless acceleration is defined as Γ = f and A are, respectively, the frequency and amplitude of the vibrations by which the system is excited, and g is the acceleration due to gravity -i.e. Γ = 1.4 corresponds to an applied acceleration 1.4 times greater than gravity. The aforementioned agreement is observed despite the fact that their simulations do not include either sliding or rolling friction coefficients, and assume a coefficient of restitution with no dependence on impact velocity. It is perhaps worth noting that earlier, similar work by Warr et al. [166] showed a significant disagreement between simulation and experiment. Interestingly, however, this is thought by the authors to be due to problems with their experimental techniques, as opposed to the limitations of their simulation model -the simulation results obtained were in fact found to agree closely with their theoretical predictions, and indeed the subsequent experimental findings of Géminard and Laroche. Studies of single-particle systems have also shown that more complex behaviour, such as period-doubling bifurcations and the transition to chaotic motion, may also be successfully reproduced by similarly rudimentary simulations; work by Tufillaro et al. [167], for instance, demonstrates how a simple model assuming frictionless, instantaneous collisions and a constant restitution coefficient can correctly predict the boundaries in frequency-amplitude parameter space, which separate the different periodic and chaotic motions exhibited by an experimental system. We next turn our discussion from single particles to the slightly more involute case of one-dimensional columns comprising multiple individual spherical particles. Once again, simulated systems assuming a frictionless environment and a constant restitution coefficient are shown to be capable of accurately emulating their experimental counterparts. The work of Luding et al. [168], for example, shows that such models can accurately recreate experimental particle distributions, scaling relations and even the presence of clustering [169], whereby repeated collisions between dissipative particles lead to localised increases in the packing density of particles, leaving other regions of the system practically empty. This faithful reproduction of experimental observations demonstrates that simulations may also faithfully reproduce dynamical processes in which particle collisions are no longer binary, but may involve multiple simultaneous particle contacts.
Simulations concerning dilute granular fluids have also been shown to accurately reproduce experimentally observed phenomena in systems with higher dimensionality. For instance, simulations of two-and three-dimensional systems are capable of reproducing the aforementioned clustering instability [169,170] and the well-known non-Gaussian velocity distributions observed in vibrofluidised granulates [171] -perhaps the two defining features of a granular gas. It is interesting to note that for the case of a highly constrained, quasi-two-dimensional system, the inclusion of frictional effects in simulations becomes highly important in producing accurate results [172], while the behaviour of fully-three-dimensional systems can still be accurately represented using frictionless walls and particles [173]. In addition to successfully reproducing the known behaviour of granular systems, simulation models have also been used to discover and characterise new phenomena unobserved in practical systems. An example of such a discovery is that of low frequency oscillations or LFOs in three-dimensional vibrated systems [174], a phenomenon first discovered using hard sphere simulations and later confirmed by experiment [175]. A similar model was also used to discover the phenomenon of sudden chain energy transfer events [176,177] in quasi-two-dimensional systems, evidence of which was latterly observed in experimental, 3D beds [178].
Having established the applicability of discrete element simulations to unary systems in one, two and three dimensions, we now relax the constraint of monodispersity, and explore the degree to which simulation models can successfully characterise the behaviour of dilute systems comprising multiple, distinct particle species. Binary beds, and indeed systems with higher degrees of polydispersity, exhibit various phenomena, such as energy non-equipartition [179,180] and segregation [16,181], which are not observed in the monodisperse case, and thus provide additional tests for the validity of discrete element models.
Feitosa and Menon [179] were among the first researchers to directly address the matter of energy non-equipartition, whereby two dissimilar particle species in the same system may possess different kinetic energies, in violation of the zeroth law of thermodynamics. Specifically, they demonstrate that binary systems of equally-sized particles, which differ in their masses and elastic properties, will exhibit decidedly unequal average energies for the two differing species. The experimental work of Feitosa and Menon [179] was later recreated in simulations by Wang et al. [182] who, unlike in many of the studies previously discussed, implement both a simple friction coefficient for particle interactions and a velocity-dependent restitution coefficient. Their setup differs from the experimental case, however, in the introduction of periodic lateral boundaries.
The simulations of Wang et al. [183] were, in the dilute regime, found to reproduce the main features of the experimental system explored by Feitosa and Menon [179], including the higher concentration of heavy particles observed near the centre of the cell, the well-mixed nature of the system on a local level, the absence of any significant clustering and, most importantly, the difference in average kinetic energy, 〈E K 〉 between particle species, with heavier particles possessing greater energy. Moreover, the experimentally observed invariance of the energy ratio between light and heavy particles with base velocity, V , and the consistency of this ratio throughout the bed's central region were also successfully reproduced in simulations.
Although the similarities in findings of the two studies discussed above strongly hint at the ability of discrete element simulations to successfully reproduce the behaviour of vibrofluidised binary granular mixtures, we nonetheless lack any direct, quantitative comparison of simulated and experimental binary systems. Indeed, to the best of the authors' knowledge, there exist few examples of such explicit comparisons in the existing literature.
We aim to address this issue -and other under-investigated systems -by including additional analysis of our own experimental and simulational data acquired from a variety of three-dimensional binary granular beds, assessing the merits and limitations of the application of DEM models to vibrated systems. Our experimental results are acquired using the Positron Emission Particle Tracking technique [34,184], while simula-tions are performed using the MercuryDEM software package [114,[185][186][187].
We look first at the situation investigated by Feitosa and Menon [179], of a binary bed of spheres with equal volumes, but differing material properties, although we focus on three-dimensional, as opposed to quasi-two-dimensional, systems. The reasoning behind this choice of dimensionality is two-fold: firstly, most granular systems in 'real world' scenarios are likely to be three-dimensional. Secondly, the dynamical behaviour of three-dimensional systems are more complex than their lower-dimension counterparts -i.e. if the behaviour of these relatively complicated systems can be successfully reproduced in simulations, then we can safely assume them to be capable of describing the simpler one-and two-dimensional cases. It is worth noting that ε may perhaps be better interpreted as an informed fitting parameter, as in reality the elasticity of a given collision is dependent on the relative velocity between colliding bodies, in addition to other, more subtle factors relating to particles' material properties and their surrounding environment . The ε values implemented are experimentally determined by measuring the pseudo-instantaneous velocity of a particle immediately before and after a colliision with a similar particle. For each particle material used, a number of such collisions spanning a range of relative velocities 0.1 v 1.0 m/s are measured and used to calculate numerous individual values of ε; the elastic coefficient for each particle species is then taken as the mean of all calculated values for the relevant species.
In fact, the only important system parameters whose experimental values are not used based on experimental data are the interparticle and particle-wall frictional coefficients, µ and µ w , respectively, which are assigned a constant value of 0.1. Interestingly, for the strongly-fluidised beds explored here, a variation in µ shows little effect on the density profiles produced, while a variation in ε exerts a significant influence, implying that collisional, as opposed to frictional, interactions are dominant in such systems. As is immediately apparent from Figure 6.3, a strong agreement between experiment and simulation is observed for the three combinations explored, strongly supporting the idea that discrete element simulations are suitable for the numerical simulation of relatively well-fluidised, binary granulates.
Not only are the observed profiles qualitatively similar for experiment and simulation, but experimentally acquired values of the systems' time-averaged centre of mass positions (and therefore mean potential energies) agree to within 10% with results obtained from simulation, a value not dissimilar to the statistical fluctuations within the experimental data used for comparison. The segregation intensity, I s -a measure to the degree of species-separation exhibited by a system [188] -is also found to agree to within 10% between experiment and simulation for the systems shown. In other words, we see that our simulations produce strong quantitative agreement with experiment, in addition to the qualitative agreement noted in our comparison of the works of Feitosa and Figure 6.3: Packing density as a function of height for a variety of vibrated, bidisperse-by-material granular beds. Data is shown for (a) a binary mixture of (light) glass and (heavy) steel particles driven with a dimensionless acceleration Γ = 6.5, (b) a binary mixture of (light) nylon and (heavy) steel particles driven at Γ = 6.5 and (c) a binary mixture of (light) aluminium and (heavy) steel particles driven with an acceleration Γ = 8. In all cases, the bed depth at rest, N L , is equal to four particle layers. Data is shown for experimental results acquired using positron emission particle tracking (blue circles) as well as MercuryDEM simulations (red diamonds) with, in each instance, open symbols representing the particle distribution for the lighter component of the bed, and solid symbols representing heavier particles. where h 0 is the bed's centre of mass height for zero driving, and h is its time-averaged, steady-state value for a continuously excited system. Data is shown for both experiment (triangles) and simulation (circles) and for two different bed depths, specifically N L = 6 (black circles and blue triangles) and N L = 3 (grey circles and red triangles). Here, gives the relative density -and, since all particles used are of equal size, mass -of the heavier (H) particle species to the lighter (L) species for the binary granular mixture corresponding to each data point. Data is shown for experimental (black circles) and simulated systems (red triangles, blue squares and green diamonds). The simulations corresponding to the triangular data points represent the 'real' combinations of materials used in experiment, with known values of inelasticity, ε, and density, ρ H,L implemented (see table 6.1). Data is shown for binary beds of steel and brass (SB), nylon and glass (NG), glass and steel (GS), glass and brass (GB), nylon and steel (NS) and nylon and brass (NB). Square and diamod-shaped data points, meanwhile, correspond to the case in which only the density ratio between particles, ρ H ρ L , is considered, with the particle elasticities held equal at values of ε H = ε L = 0.83 and ε H = ε L = 0.41, respectively. The strong agreement between experiment and simulation for the cases in which experimental ε values are utilised, compared to the considerable differences between cases in which particle elasticities (and the ratios thereof) are held constant, clearly demonstrates the importance of ε to the accurate reproduction of real-world systems using DEM models. Figure taken from our reference [189].
In fact, the agreement demonstrated in the small selection of images shown in Figure 6.3 is observed to persist across considerable ranges of parameter space. Further evidence of the ability of DEM simulations to reproduce the behaviours of experimental systems for a wide variety of system parameters and combinations thereof may be seen in our Figures 6.4, 6.5 and 6.6. Figure 6.4 demonstrates the accuracy with which simulations may predict a monodisperse system's particle distribution and, hence, centre of mass position for numerous different combinations of the driving parameters f and A as well as for differing particle numbers, N , and hence bed depths or 'layer numbers', N L . The layer number, or 'bed depth', N L is simply defined as the number of resting layers within a system (i.e. the number of layers formed when the system is exposed to zero excitation) normalised by the particle diameter, d. In Figure 6.5 we see that the agreement between simulationally and experimentally obtained values of the vertical mass centre persists also for the more complex case of binary systems, successfully accounting for variations in particle density and elasticity and, indeed, the combinations thereof for differing binary mixtures. In particular, this figure provides a stark illustration of the importance of particle inelasticity to the behaviour of granular systems, and hence the necessity of the correct implementation of this property in computational models. In Figure 6.6, meanwhile, we compare simulated and experimental segregation intensities, I s , achieved by systems in their steady states. Data is once again shown for various combinations of particle material. This image provides a fascinating insight into the importance of implementing accurate values for particles' restitution coefficients, as simulated and experimental results are observed to diverge sharply when differences in ε are not accounted for, even when all other system parameters are accurately implemented.
We have, thus far, demonstrated that DEM simulations are highly capable of accurately reproducing, both qualitatively and quantitatively, the physical properties and dynamical behaviours of relatively dilute, well-fluidised monodisperse and binary systems in one, two and three dimensions. However, as we tend toward the high-density limit in which the condition of fluidisation is no longer fully upheld and systems experience an increased prevalence of enduring contacts and a greater influence of frictional interactions, the simplified force models adopted by DEM simulations face a more stringent test.
Returning to the case of a simple, one-dimensional column -where frictional effects between particles are limited due to the restricted angle of impact between particleswe once more find that simulations can quantitatively reproduce bed behaviours; Figure  6.7, for instance, shows that -as with the dilute case -the gravitational potential energy of a system can be accurately predicted for a range of driving parameters and system sizes. Simulations are also found to accurately predict quantities such as the frequency of the periodic motion achieved by dense systems under specific driving conditions. Other studies have shown that, in addition to the global properties mentioned above, simulations may also closely reproduce the internal dynamics of a columnar system on In both panels, experimental data are represented by triangles, while squares represent simulations in which experimental density and inelasticity values are implemented and circles represent the case in which segregation is driven purely by density differences (ε H = ε L = 0.91). Filled circles correspond to experimental density ratios, once more demonstrating the considerable discrepancy between simulation and experiment when inelasticity effects are not considered. As with Figure 6 the 'microscopic' level -work by Rosato et al., for instance, has recently demonstrated that the experimentally measured [192] speed of a compression wave through an excited granular column may be correctly determined in simulation [193]; the impressive agreement observed between experiment and simulation provides a pleasing example of the ability of non-Herzian contact models to accurately reproduce physical processes even in relatively dense systems.
Studies of two-dimensional systems also show a continued strong agreement between simulations and experiment in the high-density case. Work by Yang and Hsiau, for instance, shows an impressive quantitative agreement between experimental and simulated velocity distributions [194] as well as close correspondence regarding the diffusive behaviours and granular temperatures [195] of dense, two-dimensional beds. Venturing further toward the high-density limit, where particle motion becomes less fluid-like In each panel shown, the resting gravitational potential energy, E 0 , possessed by the bed at zero energy input is subtracted from E P in order to isolate the potential energy gained by the system through the applied vibrational excitation. Figure taken from our reference [191]. and granular beds begin exhibiting crystallisation [196] and glassy properties [197], the strong similarities between experiment and simulation are still found to persist. Studies have found, for instance, that the solidus point -the filling fraction for which a granular bed 'freezes' -observed for a two-dimensional experimental system agrees quantitatively well with that pertaining to a simulated system of disks [198,199]. It is rather remarkable that such a strong, quantitative agreement is observed here in spite of several considerable differences between the experimental and simulated systems investigated, in particular the differences in particle dimensionality, with three-dimensional spheres used in experiment and two-dimensional discs modelled in simulation, and the fact that while experimental granulates are inherently non-equilibrium systems, the simulations to which they are compared correspond to an equilibrium state. Simulations are also found to successfully and accurately capture the structural configuration of dense, two dimensional systems, as demonstrated by a strong agreement between experimentally acquired and simulated radial distribution functions and shape factors [199][200][201]. In fact, the only deviations between the experimental and simulational results of the cited articles are thought to arise simply due to the imperfectly two-dimensional nature of the experimental system studied, leading to a finite number of out-of-plane collisions between particles and hence overlap in the resulting projections used to determine particle distributions. In other words, for a perfectly two-dimensional experimental set-up, one may well expect a near-perfect agreement across the full range of system parameters.
For the case of fully three-dimensional granular beds, it becomes increasingly difficult to conduct direct, quantitative comparisons of the behaviours of simulated and experimental particulate systems. Interestingly, the difficulty in comparison arises largely from the limitations of the currently available experimental techniques as opposed to the limitations of contemporary computational methods. For instance, at high packing densities, the high collision rate between particles makes it highly difficult to resolve particles' instantaneous velocities due to the limited data acquisition rates of experimental equipment [202,203], while the opacity of a densely-packed system will clearly render the real-time imaging of a three-dimensional system's interior considerably more complicated. Nonetheless, parameters such as mean squared displacement (MSD) -which reduce the noise associated with direct measurement of particle velocities -have been shown to produce a pleasing qualitative agreement between simulation and experiment in both two- [203] and three-dimensional [204] systems. Other behaviours, such as the time-dependent compaction of a granular bed and the variation of particle mobility with packing density, are also found to be qualitatively reproduced by simulations [205]. Although the existing literature -and indeed our own findings -show that DEM simulations are capable of emulating their experimental counterparts on a qualitative level, there exists, at present, little conclusive evidence to either support or disprove the ability of such simulations to provide a full, quantitatively accurate reproduction of dense, three dimensional, vibrofluidised systems.
We have thus far assessed the viability of discrete element simulation models across a broad spectrum of system parameters and dimensionalities. However, all of the research discussed up to this point pertains to granulates composed entirely of perfectly spherical particles. As we have previously mentioned, 'real' particle flows observed in industry and nature will seldom be composed of exclusively spherical constituents. However, as we have repeatedly observed from the studies discussed above, a simulation does not necessarily have to precisely match its experimental counterpart in terms of particle geometry and material properties in order to elicit an accurate, quantitative reproduction of the system's dynamical behaviours and structural properties. This fact is exemplified in our reference [208] where it is shown that the particle distributions and segregative behaviours of binary systems comprising particles differing both in density and geometry may be accurately and quantitatively reproduced using simulations of solely spherical particles. An example of the close agreement observed is depicted in Figure 6.8. The simulated density profiles shown in Figure 6.8 account for changes in particle geometry through an appropriate alteration of the diameter of the spherical particles modelled, based on the radii of gyration [209][210][211] of the non-spherical shapes used in experiment, and an alteration of the particles' restitutive coefficients, taking into consideration the additional translational kinetic energy lost to the rotational modes [212]. For the case of gaseous systems, where contacts between particles are near-instantaneous, agreement In all cases, a dimensionless bed height N L = 2.5 is used. For all images shown, particles possessing a higher radius of gyration, r g , (i.e. larger 'effective size' [206,207]) are represented by dashed lines, while solid lines correspond to lower-r g species, with black lines corresponding to experiment and yellow lines to simulation. Figure taken from our reference [208]. between simulation and experiment is strong for a variety of particle geometries and materials, across a range of driving parameters. However, for increasingly dense beds where systems become sensitive to the more specific aspects of interactions between particles as well as the differences in packing density achievable by shapes of differing geometries [213,214] this simple model is, rather unsurprisingly, found to fail.
In recent times, a number of simulation codes have been developed which are capable of directly modelling particles of various diverse geometries. Although inherently less computationally efficient than simulations involving only spherical particles, these algorithms allow a deeper investigation of the influence of particle shape on the behaviours of vibrated granulates across an increased range of parameter space.
Perhaps the simplest non-spherical geometries to numerically model are spherocylinders and ellipsoids, due at least in part to the relatively simple contact detection algorithms required [215,216]. Prior studies have shown that DEM simulations using these particular particle shapes can provide a close, quantitative agreement with experiment, with work by Pournin et al. [216], for instance, successfully reproducing the self-assembly of a bed of spherocylindrical particles experimentally observed by Villarruel et al. [217]. Even for the relatively complex case of a densely-packed, fully threedimensional system, a comparison of the two studies shows quantitative agreement in terms of the initial and final packing densities achieved and even the timescales corresponding to the ordering processes within the systems explored. Beyond the constraints of ellipsoidal/spherocylindrical geometries, there exists an apparent lack of DEM-based research focussing on vibrated systems with more complex non-spherical constituents, despite the existence of numerous such studies pertaining to sheared beds [218,219] and gravity flows [183,220], for instance. Recently, however, the multi-sphere method [221], which has proved highly effective in the modelling of other system geometries [222][223][224], has begun to be applied to vibrated and vibrofluidised beds. Using the multisphere approach, each non-spherical particle within a system is composed of a number of smaller objects, creating 'composite particles', the modelling of which is often more computationally efficient than for complex shapes represented using single particle models [225,226]. Work by Chung et al. [227], for example, demonstrates, for the case of simple paired particles, an excellent, quantitative agreement between experiment and simulation, providing strong support for the validity and applicability of the multi-sphere model in vibrated beds. Pei et al. [228], meanwhile, implement more complex elongated particles comprising spheres of differing sizes. They study the effect of geometry on the transfer of electrostatic charge between particles, finding their results to agree qualitatively with experimental expectations.
Qualitative predictions arising from numerical studies using the Monte Carlo method [209] to simulate the segregative behaviours of particles possessing still more complex geometries also agree well with experimental observations [208]. Although these simulation methods can at best be described as semi-quantitative, the observed correspondence between experiment and simulation nonetheless shows promise.
In summary, we have demonstrated that numerical simulations of vibrated and vibrofluidised granular beds can provide quantitatively accurate models of the behaviour of these systems over a wide range of parameter space. The models used are capable of recreating numerous complex behaviours exhibited by particulate systems spanning one, two or three spatial dimensions, existing in solid-like, liquid-like or gaseous states and comprising particles possessing various diverse material properties. While there remain certain cases where the applicability of particle simulations is yet to be fully verified, with the continuing improvements in processing power, in addition to the further refinement and development of the computational methods and contact models associated with DEM simulations, it is likely that these remaining issues will be resolved in the not-too-distant future.

Flows over inclined channels (chute flows)
Similar to vibrated systems, inclined channels (chutes) are also commonly employed in various industrial processes, for example in the mining, pharmaceutical and foodprocessing industries, which handle particulate media [2]. Besides their numerous industrial applications, a range of geophysical events (including avalanches, landslides, debris flows, see [229][230][231]), can be emulated utilising this simple flow geometry. Hence, a thorough understanding of the relevant granular dynamics is not only vital for designing efficient handling equipments, but also for developing accurate continuum formulations used to predict these hazardous natural events 5 . Attributable to the inherent nature of granular materials to exhibit a range of diverse behaviour, they are broadly categorised into three regimes: solid, liquid, and a gaseous (see Fig. 6.9). Vibrated systems (section 6.3.1) can, under certain circumstances, fall into the category of gases because the particulate media is dilute and the grains interact through binary collisions. However, depending on the level of energy input and the chosen dissipation parameters, these systems can appear to be solid-or liquid-like as well. In contrast to vibrated systems, granular flows over inclined channels fall in the category of intermediate liquid regimes, where the flows are rather rapid and dense -high volume fraction -and the momentum exchange and energy dissipation is due to both collisional (instantaneous) and frictional (enduring) contacts between the particles. This section focusses solely on these dense liquid flow regimes. Essentially, a granular material is a conglomerate of multiple discrete, solid constituents, which are typically dispersed in a medium. In rapid dense flows, the particles are much denser than the interstitial fluid with low dispersion (closely packed). As a result, the interstitial fluid, e.g. air, is often neglected as it has negligible effect on the bulk flow behaviour. Despite such simplifications, the behaviour of these dense rapidly flowing dry particulate media is yet to be completely understood.
In order to gain an improved understanding, studies have utilised a varied range of experimental setups -plane shear, annular shear cells, vertical chutes, rotating drums, heap flows, and inclined channels (our current interest) -in which the granular material is subject to simple or planar shear. For the schematics of these setups, see MiDia [233]. By employing the available state of the art experimental techniques, such as particle imaging velocimetry (PIV) [115], particle tracking velocimetry (PTV) [234], force sensors [235,236], acoustic probes [237], tracked transmitters [238], capacitance probes [239], optical imagery [240], digital imaging [241][242][243], refractive index methods (RIM) [244], X-ray tomography [245], magnetic resonance imagery [246], studies have not only determined the kinematic properties, but also the rheological behaviour of dense granular flows. However, in 3D, no contemporary experimental technique has the ability to directly determine the local or global stresses generated in such flows. In contrast, discrete element simulations are able to provide all the useful information necessary to fully understand the kinematics and dynamics of a system -hence this indispensable need for cooperation and agreement between experiments and (discrete element) simulations.
Although simple, understanding dense inertial flows has been and still is a daunting task. We begin with a relatively simple case of monodisperse flows over both smooth and rough inclined channels. Interests in understanding these flows date back to the 70s' where several studies [247][248][249][250][251] focussed on determining the rheological properties (constitutive laws) of mixtures constituting sand and glass beads. Later Hwang and Hogg [252] and Brennen et al. [253], focussed on examining the diffusive remixing in these flows and illustrated the possibility of granular jumps or bores, respectively. However, due to the lack of intrusive experimental techniques to precisely determine the flow velocity, solid volume fraction and granular temperature, many studies utilised computer simulations as an alternative approach, see [254,255] and references therein. Campbell et al. [254,256] were one of the first to compare their particle simulations of twodimensional unidirectional flow of inelastic circular cylinders with experiments. They used a hard-sphere DEM approach and observed their results to qualitatively match the findings of Augenstein & Hogg [250] and Bailard [257], but not of Savage [251] and Ridgway & Rupp [247]. While the majority of the studies focussed on understanding the fundamentals [258][259][260][261][262][263][264] and developing continuum formulations [145,265,266], it was in the early 90s' where discrete element simulations were employed to analyse the two-dimensional [267][268][269] and three-dimensional [270] flow behaviour of inelastic frictional particles. However, most of these simulations either contained very few particles -compared to todays standard - [267][268][269] or were not fully steady [267].
In the early 2000s, a few of the calibration steps illustrated in Fig. 6.2a were pursued to investigate flows of glass beads over a rough inclined channel [271]. Using a combination of experiments and discrete element simulations, Hanes and Walton [271] studied the effects of basal and side-walls roughness on the flow structure and its dynamics. In their experiments, steady flows were observed for a range of inclinations. Furthermore, for two inclinations, the particle velocities (both near the side-walls and at the free-surface), the mass flow rate and the flow depth were found to be in relatively good agreement with discrete element simulations. Furthermore, their DEM simulations utilised the contact parameters for glass beads which were determined using single particle experimental measurements carried out by Lorenz et al. [272]. Similarly, several other studies [35,149,273,274] have used DEM to simulate both steady and unsteady flows of cohesionless particles (glass beads) and found them to be in good agreement with experimental observations. Silbert et al. [273,275] investigated the rheology of fully three-dimensional, steady, dense flows of cohesionless frictional spheres over rough inclined planes. The channels were made rough by sticking particles on the channel surface, over which the material flows. Among the areas studied were the effects of parameters, such as the inter-particle friction coefficient, inelasticity, flow heights and channel inclinations. Steady flows were found to exist over a wide range of inclination angles and heights confirming the experimental observations of Pouliquen [148]. For different flow heights (thick flows), Silbert et al. [273] found the flow properties, such as velocity and strain-rate, to be sensitive to the interparticle friction coefficient. Additionally, for thick enough flows the velocity profiles and rheology were shown to obey Bagnold scaling. However, when gradually reducing the flow height, Silbert et al. [275] observed transitions from Bagnold rheology to linear velocity profiles, consistent with the experimental observations of Lemieux & Durian [276] (avalanches), Ancey [277] (linear velocity profiles) and Pouliquen [148] (Bagnold scaling). Furthermore, Silbert et al. [275] observed that both thin and thick flows obey a simple scaling law consistent with the experiments of Pouliquen [148]. In addition to the range of parameters -such as the flow height, angle of inclination, etc. -considered to understand their effects on these flows, Weinhart et al. [35] and Thornton et al. [149] also systematically investigated the effects of varying basal properties. By changing the basal particle diameter or inter-particle friction between the flowing particles and fixed base particles, [35,149], they studied the effects on the effective macroscopic friction coefficient. For varying basal particle diameters, Weinhart et al. [35] found their simulation results to be in good agreement with the empirical law determined by the experiments of Pouliquen & Forterre [148,278]. Not only did Weinhart et al. [35] obtain a DEM validated empirical law for the effective macroscopic friction coefficient, but they also determined other closure relations, such as for ratio of the stress in the downslope and normal direction (flow depth), K , mean flow density,ρ and the shape of the velocity profile, α, -essential for solving the continuum shallow-layer granular continuum model, e.g. as shown in Chapter 2 or references [145,147]. Similarly, for varying basal contact friction coefficients, Thornton et al. [149] presented a modified empirical law [278], which could be used in continuum models to simulate flows over varying base types, i.e. rough or smooth. In addition, by choosing the right inter-particle collision model and contact parameters, recent experimental studies have further validated discrete element models of monodisperse glass spheres over nonrotating [279] and rotating smooth inclined channels [116,280], thereby illustrating once more the current ability and future potential of discrete element simulations to emulate and help understand monodisperse flows over both smooth and rough inclined channels. For more examples and a detailed review concerning other applications of discrete element simulations, see [74,233].
As stated earlier, real-world flows comprise particles of varied sizes, densities, shapes and many other characteristic features. Attributable to these varied features, particles in dense rapid flows (also in vibrated systems) tend to arrange themselves in distinct spatial patterns [10,[283][284][285][286], which is termed as particle segregation. Despite many differences in particle attributes, differences in size and density are considered to be the most important, see [285,286] and references therein. Similar to the interest shown in understanding monodisperse flows, several experimental studies have attempted to address the issue of size-based particle segregation in rapid dense flows over inclined channels [265,279,[287][288][289][290][291][292][293]. It was in the late 90s' when a two-dimensional discrete element simulation of inelastic disks was carried out to quantitatively describe particle segregation [294]. However, no quantitative agreement was achieved with the experiments of Savage & Lun [265] because of the different mixture compositions considered. While many studies, e.g. see Tunuguntla et al. [281], Fan and Hill [295], Schlick et al. [296], focussed on developing accurately predicting continuum segregation models using DEM, Thornton et al. [186] illustrated the need of a frictional tangential force collision model for producing steady flows with strong segregation as seen in experiments. Interestingly, similar results were also produced by simply adding tangential dissipation to the contact model [186]. Although they were able to observe realistic flow behaviour, further investigation is necessary to understand the effects of the chosen contact models on segregation. Recently, by a systematic fine tuning (parameter calibration) of the contact model, Bhattacharya & McCarthy [279] showed quantitative agreement with experiments. While purely size-based segregation has been a topic of interest for several years, Tunuguntla et al. [281] (Chapter 4) investigated both size and density effects on Figure 6.11: The above phase-plot is constructed using several fully three-dimensional periodic-box DEM simulations of bidisperse mixtures, comprising of two types of particles. Using a simple force model, homogeneously mixed bidisperse mixtures, for a range of size-and density-ratios, are allowed to evolve till steady state is reached. Given s = s 2 /s 1 and ρ = ρ 2 /ρ 1 as the size-and density-ratios, the above illustration plots the normalised relative vertical center of mass of type-2 particles, D com = (z 2 com − z b com )/z b com , with z 2 com the vertical center of mass of type-2 particles, and z b com the bulk vertical center of mass. As an example, for a purely size-based mixture ( s, ρ) = (1.4,1.0), we have D com = 0.5, implying that the large particles rise towards the free-surface whereas the small particles settle near the base of the flow. Additionally, the solid line denotes the line of weak segregation predicted from the continuum theory model of Tunuguntla et al. [281] whereas the dotted-dashed line is the weak segregation line predicted analytically by Jenkins & Yoon [282] using kinetic theory. For more details see Tunuguntla et al. [281].
segregation, see Fig. 6.11. Although no direct quantitative comparisons were made with experiments, qualitative agreements were found with the experiments of Felix & Thomas [288]. These results imply that, wih proper calibration of DEMs, one should be able to effectively emulate the experimentally observed particle phenomena on a quantitative level.
Although studies have shown particle shape to have significant effects on particle dynamics, there exist very few studies concerning the experimental validation of DEMs of non-spherical particles flowing over inclined channels. Nevertheless, by proper calibration, Vu et al. [297] illustrated the ability of DEMs to simulate flowing ellipsoids in an attempt to simulate the dynamics of flowing soybeans, their results being further validated by experimental findings. However, additional validation of particle simulations is required to comfortably apply them to the analysis of the complex flows of non-spherical particles.

Summary and Conclusions
Through analysis of the extensive literature related to the computational modelling of granular and particulate flows, we provided an overview of the extent to which current discrete element simulations may successfully emulate the behaviour of dynamic experimental systems, focussing specifically on vibrofluidised beds and chute flows -two systems with direct relevance to numerous natural and industrial processes.
We carefully examined the specific circumstances under which computer models may -if correctly implemented -be expected to accurately recreate experimental systems for two commonly explored types of granular flow. We show that simulations of vibrated systems using the discrete element method are capable of quantitatively predicting the behaviour of granulates in one-, two-and three-dimensional systems, for the case of monodisperse beds and bidisperse mixtures, and for systems in dilute gaseous states, liquid-like states and densely-packed solid-like states. Discrete element models of chute flows, meanwhile, are found to show strong quantitative agreement with experiment for monodisperse flows across flat, bumpy, frictional and smooth surfaces for all dimensionalities. The segregative behaviour of binary systems, which are bidisperseby-size, have also shown to be correctly simulated using the discrete element method, although agreement between simulated and experimental systems comprising particles of differing densities remains to be directly tested.
Perhaps the most pressing open question for both vibrated granular systems and chute flows is whether the effects of particle shape can be adequately captured by simulation models. At present, although quantitative agreement has been observed for very simple systems, and qualitative trends reproduced for more complex cases, there exists little direct evidence to concretely confirm or deny the suitability of current models to the task of modelling non-spherical systems. Considering that the majority of particles involved in real industrial and natural processes possess non-spherical geometries, this is most certainly an issue worthy of direct, systematic research in the future.
In conclusion, the collected works discussed here demonstrate the ability of discrete element simulations to quantitatively emulate experimentally observed phenomena exhibited by dry granular media in both vibrated beds and liquid-like regimes (chute flows). In doing so, we show that such simulation models provide a potentially power-ful tool with which we may better understand and predict the behaviour of granulates in these commonly used flow geometries. The successful application of DEM models to these canonical systems suggests that they may indeed prove equally valuable in other granular systems.
So many worlds, so much to do, so little done, such things to be.

-Alfred, Lord Tennyson
The focus of this thesis concerned modelling the dynamics of rapid dense granular materials flowing over inclined channels, using in-depth theoretical analysis, discrete particle simulations (DPMs) and an accurate micro-macro mapping technique. By a thoughtful combination of each of these individual elements, a beautiful blend has been established among different scales, i.e. from particle to continuum. Overall, a sincere effort has been made towards developing this blend, which is able to help understand and emulate these phenomena-rich inclined channel flows.
As a primary step in our investigation, we begin in Chapter 2 by considering shallow mono -disperse flows over a channel with converging sidewalls. By shallow, we implied that the ratio of the characteristic length scales in the normal (H) to the streamwise direction (L) is small, H/L << 1. Using the same shallowness argument in the cross-slope channel direction, we derived a novel one -dimensional (width-and depth-averaged) granular hydraulic theory from the depth-averaged shallow granular model. For closure, the model utilised an empirically determined DPMs validated constitutive law. More importantly, the closure law accounts for the existence of steady uniform flows for a range of channel inclinations. On solving the one-dimensional model for steady flows, the theory predicted the existence of multiple flow regimes. The flow regimes ranged from regular steady smooth flows to the ones with steady jumps in their height or velocity profiles. To strengthen the value of this model, Chapter 2 further illustrated a much required verification step, where we numerically solved the two-dimensional shallow granular model using a continuum solver, discontinuous Galerkin finite element method. For supercritical flows, we compared the flow height obtained from both the one-and twodimensional theories. Surprisingly, the profiles matched well, despite the presence of 145 oblique granular jumps across the contracting channel. However, further verification is to be done for, both, flows with jumps and subcritical flow regimes.
The majority of inclined channel granular flows involve polydisperse mixture rather than an ideal monodisperse material. Thereby, constructing a continuum description to describe these flows is not only enticing, but also a formidable challenge. Several studies have utilised DPMs as an alternative to experiments to understand the particle scale dynamics. However, an accurate mapping of the particle scale mechanics onto a continuum macroscopic field is still an area of ongoing research.
Chapter 3 comprehensively presented a generic framework for an efficient and accurate micro-macro mapping technique, called coarse-graining, for polydisperse mixtures comprising of convex shaped constituents, e.g. spherical particles. More importantly, the method presented is valid for any discrete data, e.g. particle simulations, molecular dynamics, experimental data, etc. However, for the purpose of illustration, Chapter 3 considered data generated from discrete particle simulations of bidisperse granular mixtures flowing over rough inclined channels. In order to obtain continuum macroscopic fields, for any stationary or transient particulate system, it is essential to choose a proper spatial coarse-graining scale, irrespective of the chosen coarse-graining function. Thus bringing us to the question what do we mean by a continuum description? Chapter 3 comprehensively answers this question, through the example of bidisperse mixtures flowing over inclined channels. We showed that when one chooses a length scale smaller than the continuum length scale, the resulting coarse grained data will still show individual particles; these are not continuum fields. On the other hand, if one chooses a large coarse-graining scale, it will smear out the macroscopic gradients. Thence, the results will be strongly dependent on the chosen coarse-graining scale. However, Chapter 3 showed an existence of a range of spatial and temporal coarse-graining scales in which the continuum fields obtained are independent of the chosen coarse-graining scale. It is this length and time scale that must be utilised for an efficient micro-macro transition. Once an optimal scale to coarse grain is determined, more importantly, Chapter 3 illustrated the utility of coarse-graining the discrete particle data for both steady and time evolving particulate mixtures. With this availability of accurately constructed continuum fields, now one could methodically develop and calibrate an efficient continuum model.
Albeit the continuum description in Chapter 2 attempts to model the dynamics of monodisperse flows over inclined channels, it is still unable to capture other complex phenomena such as particle segregation. Before employing the efficient mapping technique, in Chapter 4 we formulated a mixture theory based segregation model for bidisperse mixtures varying in, both, size and density. The model is based on the current understanding of the bidisperse segregation dynamics, which is based upon an idea/assumption reflecting on how the total pressure/bed load during the flow is distributed among the two mixture constituents, i.e. the so called pressure scaling functions. The developed formulation is built upon an existing size-segregation model, and is applicable to both shallow (linear velocity profile) and thick (Bagnold profile) flows. For linear velocity profiles, the current model has also been analytically solved and compared with the few of the existing size-segregation models. Besides predicting the extent of segregation for range of size-and density-ratios, more importantly, the theory also predicted zero or weak segre-gation for a range of size and density ratios which was further confirmed using DPMs of spherically shaped particles.
Chapter 4 considered an existing form of pressure scalings. In Chapter 5, we investigate this in more details using the coarse-graining technique of Chapter 3. For simplicity, we considered a purely size-based segregation model, which was built upon a simple pressure scaling function containing an unknown parameter. Not only did we determine this unknown material parameter using coarse-graining technique but, more importantly, we also found out that the complete size-and density-based segregation in any flowing particulate mixture is an effect of the generated kinetic stress rather than the contact stress.
Chapter 3, 4 and 5, illustrated how one could blend continuum models with DPMs, to generate predictive continuum models. However, an important point that has not been addressed so far, is to see if the DPMs can actually emulate reality. As a consequence, in Chapter 6 we investigate, how well particle simulations capture the reality.

Outlook
Given the developments illustrated in this thesis, we present here some prospects for future work: ♣ The novel one-dimensional (1D) and two-dimensional (2D) shallow granular model predicts the flow dynamics for monodisperse flows alone. To enhance the strength of the developed one-dimensional model, a proper validation -either through experiments or DPMs -would definitely benefit the theory.
♣ Additional value can be added to the 1D and 2D shallow granular model, if one does properly couple the two-dimensional shallow granular model (SGM) with the bidisperse segregation model. The 2D SGM constructs the velocity profile and provides it as an input for the segregation model. However, a suitable closure law should also be determined to take the bidisperse nature of the mixtures into account. Although, a recent study has shown the Pouliquen friction law to be valid for bidisperse mixtures varying in size alone, further validation needs to be done to check if the closure law is valid for bidisperse mixtures varying in, both, size and density. The developed model can then predict phenomena such as granular fingering, more accurately.
♣ As a step towards developing calibration and validation tools, i.e. from particle to continuum scale, the micro-macro transition using the coarse-graining technique needs to be first applied to analyse polydisperse mixture dynamics. In addition to this, it further needs to be extended to non-convex shaped particles. Suitable changes need to be accounted for incorporating non-convexity into our existing theory which already considers polydispersity without any loss of generality.
♣ Additionally, the coarse-graining tool needs to be generalised, i.e. currently our method carries out the spatial coarse-graining where one could choose different coarse-graining functions. These are Heaviside, Gaussian and Lucy functions. However, this is not true with respect to temporal coarse-graining. At the moment, we assume Heaviside function to be our temporal coarse-graining function. This feature definitely needs to be extended to consider other variations.
♣ The prediction of zero segregation from the bidisperse segregation model was confirmed using bidisperse mixtures with each constituent' solid volume fraction as 50%.
Further study needs to be carried out for checking the effects of solid volume fraction on the zero-segregation prediction.
♣ To enhance the predictability of the bidisperse mixture theory models, one definitely needs to determine the correct form of pressure scalings. This still is an area that is vastly unexplored. The second suitable extension to the bidisperse segregation model would be to include diffusive remixing correctly. Several continuum descriptions, modelling segregation, assume a certain way of capturing diffusive effects. However, detailed study needs to be carried out in order to determine a correct continuum description for diffusive remixing.  Given we have a shallow flow, i.e. the ratio of the characteristic length scales in the normal (H ) to the streamwise direction (L) is small, H /L << 1. The two-dimensional depth-averaged model for shallow granular flow is derived using the Cauchy's equations for balance of mass and momentum, where ρ is the material density, u = [u, v, w] T is the bulk velocity, σ is the stress tensor and b is the body force, e.g. gravitational force. The above balance laws are to be satisfied by all solids, liquids and gasses. In addition, D/D t = ∂/∂t + u ·∇ is the material derivative operator.
In order to arrive at the shallow granular model, a series of assumptions or approximations are made. The material is assumed to be incompressible with a constant material density ρ 0 and gravity being the only external body force. As a result, the simplified Cauchy's equations are restated as where g is the acceleration due to gravity and θ is the inclination of the coordinate system with the horizontal, see Substituting suitable scalings and boundary conditions, see [1,2], the above Cauchy's balance equations are depth-averaged to result in the two-dimensional dimensionless shallow granular equations (A.4) where u and v is the depth-averaged velocity component in the down-and cross-slope flow direction, α 1 and α 2 denote shape factors, i.e. iff denotes a depth-averaged quantity, thenf 2 = (f ) 2 . Thereby, we haveū 2 = α 1 (ū) 2 andv 2 = α 2 (v) 2 . For approximately uniform velocity profiles, α 1,2 ≈ 1. In addition, K represents the material constant denoting the stress anisotropy, i.e. σ xx = K σ zz . In the simplest case K = 1. Furthermore, ǫ is the aspect ratio H /L and µ is a dimensionless friction coefficient representing frictional effects in an average sense. For more details concerning the derivation of a shallow granular model, see [1,2].

A.1. Deriving a one-dimensional shallow granular model
For deriving the 1D shallow granular model, we restate the depth-averaged granular model in its dimensional form. Assuming constant/no basal topography (b(x, y, t )=const. or 0) and shape factors α 1,2 = 1, we have The subscripts t , x and y denote partial derivatives. Applying the approximate 1D hydraulic analysis, we shall average the flow quantities across the cross section of a channel gradually varying in width. The channel sidewall is symmetric with respect to the x-direction and remains constant in time, see right sidewall are defined as y = ± W (x,t ) 2 . The kinematic boundary conditions must hold at the left and right sidewall. Let the left sidewall be denoted by superscript 'l', it implies Similarly for the right sidewall, denoted by superscript 'r ' Assuming that the boundaries are impermeable, for fixed boundaries we obtain a more familiar boundary condition u · n = 0, where n is the outward normal at the fixed boundary.

Mass Balance
Applying Leibnitz' rule on each term of (A. 8 Using the boundary conditions, (A. 6) & (A.7), in the width-averaging of the third term results in The width-averaged valuef of a variable f is given byf = (1/W ) Adding all the above results, dropping the bars, we obtain a width-averaged mass balance equation .20) Momentum Balance

Mass Balance
After dropping the apostrophes, the non-dimensional 1D shallow granular equations are From above relations, (A. 24) A. 1 Substituting h = QF l W F 2/3 , results in 25) A. 1

.5. Regularisation
The limit problem is stated as where x cnew is the new contraction exit determined during regularisation. Below we solve the limit problem for both inviscid and viscid cases.

Inviscid case
The slope of the contractionγ and the normal to the contractionˆ n are given bŷ Figure A.3: x c in the initial x-location of the contraction exit. x r is the x-location for the center of the circle used to regularise. For inviscid case, after regularisation x r turns out to be the new x-location for the contraction exit.
The equation of a line with position vector x c W c and direction vectorˆ n is As x r ,W r lie on the line (A.28) The value of λ r is determined by considering the distance between (x c ,W c ) and (x r ,W r ), which is using (A.30), gives us With λ r determined, the center of the circle for regularisation purpose is computed using (A.30). By using the equation of a circle, we have

Viscid case
Eqn. (A.25) can be restated as Using (A.38), To compute the Taylor series expansion of (A.41), we use the following expansions (A.42) By utilising the above expansions, we compute the Taylor expansion for the term below as (A.43) Thus, using A.43, a Taylor expansion for the term below is computed as follows v (F l Q) 2/3 1  47) Substituting all the above approximate expansions, the Taylor expansion for   A mixture theory segregation model is formulated using a postulate which states that the constituents of the mixture simultaneously occupy the same point in space and time. This leads to overlapping fields like partial pressure p ν , partial density ρ ν , partial velocity u ν , etc., corresponding to each of the constituent ν = 1, 2. For a mixture comprising of several constituents, each of the constituent satisfies the fundamental balance laws of mass and momentum. These laws are stated in terms of partial fields, as below where g is the acceleration due to gravity vector, g = g t , 0, −g n with g t = g sin θ and g n = g cos θ. The drag force β ν is the net effect of the tractions across the interface of The contents in this appendix is to be considered as a supporting material for Chapter 4. two constituents. As the species percolate or lever during kinetic sieving (segregation mechanism), the species experience interspecies friction. Experiments showing the kinetic sieving phenomenon suggest that the process is analogous to percolation of fluids through porous solids. Therefore, the interaction drag or interspecies friction is assumed to have the form where c is the coefficient of the interspecies drag,γ is the shear rate and φ ν is the volume fraction that denotes the percentage of the local volume occupied by each constituent.
In reality the mixture flow consists of three-constituents, the additional component being the void space between the particles. However, the effect of including this space is minor and is discussed in more detail in Thornton et al. [1]. The first term of (C.2) is the interspecies surface interaction force, which when incorporated within the momentum balance law leaves behind φ ν f ν ∇p. This enforces the percolation process to be driven by intrinsic pressure gradients. The second term defines the linear rate-dependent resistance to relative motion. With the definitions of partial pressure and drag force at hand, we have ∇p ν = ∇(pφ ν f ν ) = p∇(φ ν f ν ) + φ ν f ν ∇p. The momentum balance law for each individual species (C .1) 2 , using the mass balance law (C .1) 1 , is then restated as

C.2. Boundary conditions
By dropping tildes for simplicity, we have Basal boundary condition Using above scalings (C.4), i.e. non-dimensionalisation leads to On dropping tildes for simplicity, we have

Boundary asymptotic
Besides carrying out the asymptotic analysis on the mass and momentum balance equations, the asymptotic expansions are substituted into the expressions valid at the boundaries as well.
Using the asymptotic expansions, we have (C.38)

C.4. Segregation governing equation
Defining φ 1 0 = φ, the continuity equation is stated below, Utilising the results obtained from the above asymptotic analysis, where the partial downand cross-slope velocity is, in the leading order of the aspect ratio, equal to the bulk down-and cross-slope velocities u 1 0 = u 0 , v 1 0 = v 0 , then (C.26) is restated as  For a purely density driven segregation, s 1 = s 2 , the governing equation is

C.5. Summary
From the x-and y-momentum equations, O(1/δ 2 ) implies that the partial velocity in the down and cross-slope is equal to the bulk down and cross-slope velocity, 47) their patient perusal of this thesis and, more importantly, for the support and belief when needed. I would also like to thank Marielle, Linda and Sylvia, secretaries of the groups Mathematics of Computational Science (MACS) and Multi-Scale Mechanics (MSM) under MESA+ institute, for taking care of all the administrative work related to the project and my employment at the university. Without them, things would definitely not have been simple. Besides my heartfelt appreciation for the supervision and administration, the decent quality of research -presented in this thesis -would not have been possible without a motivating working environment and a healthy social life. I thank the members of MSM and MACS groups for providing me with the former. Regarding the later, life in this small town (Enschede) would definitely have been boring without a lively ambience. May it be partying or chilling (as they call it these days), this small university town has certainly made my stay here joyful. And of course, all this would have been a mere fantasy without the tons of people I have met here. Naming each one of them would be a daunting neverending task and hence skipped.
Finally, no story is complete without the umpteen support and love from the family. Thank you amma, nana and Pooja for believing in me all these years, and bringing along some real food during your visits here. Also special thanks to the van den Toorn's for being the family away from the family. And last but definitely not the least, although words are not sufficient, thank you Merel for being there and bearing all of my drama.
Thank you all for making Enschede a lively and loving place.

Summary
This thesis focuses on modelling the dynamics of dense granular materials flowing over an inclined channel, utilising a continuum description. In the process of understanding and developing this, besides continuum modelling, the thesis also exhibits the utility of discrete particle simulations (DPMs), and the need for developing an accurate micromacro mapping technique.
As most of these dense gravity-driven flows are shallow in nature, for monodisperse mixtures, Chapter 2 illustrates the formulation of a novel one -dimensional (width-and depth-averaged) shallow granular model. Using this model, we not only predict the flow dynamics -flow height, velocity and granular jumps or shocks -but also shows that one can forecast the existence of multiple steady states for a given a range of channel inclinations. However, in reality, the majority of flowing particulate mixtures are known to comprise of particles with varied physical attributes, i.e. they are bidisperse 1 or polydisperse 2 . Thereby, as a step towards understanding the associated flow dynamics, and developing improved continuum models, several studies presented in this thesis have utilised discrete particle method. DPMs provide a plethora of information at a particle scale, such as particle position, velocity, interaction forces or stresses. In order to accurately map the particle scale mechanics onto a macroscopic continuum scale, Chapter 3 comprehensively presents a generic framework for an efficient and accurate micromacro mapping technique for polydisperse mixtures comprising of convex shaped particles, e.g. spheres. More importantly, the method presented is valid for any discrete data, e.g. particle simulations, molecular dynamics and experimental data, for both steady and unsteady configurations.
Before employing the efficient mapping technique of Chapter 3 to its full capacity, based on the current understanding of bidisperse segregation dynamics, we formulate in Chapter 4 a mixture theory segregation model for bidisperse mixtures varying both in size and density. The developed formulation is an extension to an already existing size-segregation model, and is applicable to both shallow (linear velocity profile) and thick (Bagnold profile) flows. Besides predicting the extent of segregation, the theory also predicts zero or weak segregation for a range of size and density ratios, which was further benchmarked using DPMs. Although, we developed an efficient continuum sizeand density-segregation model, a detailed study is to be implemented in order to determine more accurate pressure scalings and further extend it to polydisperse mixtures. As a stepping stone, towards determining these pressure scalings, in Chapter 5 we give an example application of the micro-macro mapping technique (illustrated in Chapter 3). For simplicity, we consider a purely size-based segregation model, which was built upon a pressure scaling function containing an unknown parameter. Not only did we determine this unknown material parameter but, more importantly, we also found out that Summary the complete size-and density-based segregation in any flowing particulate mixture is an effect of the generated kinetic stress, rather than the contact stress. The current form of the existing scaling functions is, however, still an active area of research, which definitely needs further attention and care.
Chapters 3, 4 and 5, show how one can mix and match continuum models with DPMs using an efficient coarse-graining method. However, it is still vital to see if the DPMs can actually emulate reality. As a consequence, we illustrate in Chapter 6, how DPMs can be used as a suitable alternative to experiments using two commonly used DPM experiments.