Size segregation in compressible granular shear flows of binary particle systems

This paper deals with the modelling and simulation of segregation in granular materials. The basis is a hydrodynamic model for granular material flows, which is extended to capture the dynamic process of segregation in shear flows of systems with small and large particles. The granular flow equations consist of a set of compressible Navier–Stokes-like equations as well as an equation for the granular temperature. With the help of the granular temperature equation, the granular flow equations are able to cover a wide range of regimes, starting from dilute to arresting flows. However, this paper focuses on dry granular shear flows. It extends this hydrodynamic system in a dense shear flow regime by a segregation equation using the framework of mixture theory. Special focus is lain on the segregation direction. A procedure from mechanics is adapted to obtain the segregation direction from the granular flow system independent of the choice of the coordinate system. In particular, this is done in three-dimensional space. Due to the compressibility of the granular flow system and the structure of the derived segregation equation, solving the segregation equation requires special numerical treatment. Therefore, a suitable numerical scheme is presented which prevents the system from reaching unphysical states.


Introduction
Granular materials consist of macroscopic particles of different size and kind. Examples range from avalanches and dunes in geophysics, grain filled silos in agriculture to powders in the cosmetic and pharmaceutical industry. They are omnipresent in humans daily lives. Many industrial and chemical processes have to deal with granular materials. The worldwide annual production of grains and aggregates reaches 10 billion metric tons and their processing consumes 10% of all energy produced on earth. It is the most manufactured material in industry after water [7]. For these reasons granular materials have been the subject of intensive engineering research for many decades, but fully understanding their dynamic behaviour still poses a major challenge to engineering science and physics. The natural and industrial processes described above can be quite complex, since granular materials can behave similar to either fluids or solids. For the aforementioned processes, the different regimes occur at the same time or consecutively. Due to the macroscopic size of the discrete particles, Brownian motion has no relevance for the collective behaviour. The dissipative forces acting on the particles, such as inelastic collisions and friction, lead to different material properties than conventional fluids, solids or gases. For a granular system to remain active, it needs to gain energy from external forces (gravity, electric or magnetic fields), shear or vibration. Thermal fluctuations are insufficient to move grains and therefore do not play a role for granular systems. Consequently, granular materials exhibit metastable steady states far from equilibrium [1]. All these aspects contribute to the difficulties one encounters in the field of modelling granular flows.
To simulate the behaviour of granular materials, for researchers, the field of granular physics is still a mixture of different modelling tools, concepts, and phenomenological theories. After all, there is no general hydrodynamic theory as the Navier-Stokes equations (NSE) for simple fluids, which can model a granular system with the same accuracy. Quite good results can be acquired from discrete models like the discrete element method (DEM), but they are too costly to simulate entire industrial processes. With kinetic theory there exists a unified description for dilute systems of rapid grains, but it is not applicable in general for granular systems approaching close-packing density. However, it has been shown that extensions to kinetic theory can be derived that are applicable to denser flows [22].
In case of granular systems composed of different grain sizes, segregation effects may occur during flow. The segregation effect might be very helpful in the mining industry for mineral processing, but it leads to problems in many other industrial areas. In the pharmaceutical and the food industry, mixing processes using rotating tumblers are very common. A uniform mixing is also desired in the production of ceramics or plastics. Here, segregation can degrade the quality of the products [40].
Since size segregation is the most dominant mechanism compared to segregation effects due to particle density or particle shape, the focus lies on the topic of size segregation. In recent years, the starting point for the derivation of most segregation models is to look at a simple process. Typically, granular avalanches of small and large particles flowing down an inclined plane serve as a common example. For the mathematical description of the segregation process, an advection-diffusion equation for the solid volume fraction of the small particle phase ∈ [0, 1] with a certain non-linear structure of the flux term is established in literature. In the most publications, it can be given by where u defines the bulk velocity in x-direction, q the segregation velocity, and D the diffusion coefficient. A first one-dimensional model, which ignores the transport with the granular material and a flux term with = 2 , was given by Bridgwater et al. [3]. Later models are similar to equation (1) but with a flux term where = 1 , like in the work of Dolgunin and Ukolov [6]. Since 2005, mixture theory, which is a theory to model multiphase systems using the principles of continuum mechanics, has often been used to derive segregation equations. A first rigorous derivation of an equation like (1) was given by Gray and Thornton [19]. The authors concentrated on the advective segregation process without diffusion. An extension including the diffusive remixing term was introduced by Gray and Chugunov [18]. Further extensions were introduced later on concerning multi-component segregation by Gray and Ancey [16], asymmetric flux modelling by Gajjar and Gray [10], or combined particle-size and particle-density segregation by Tunuguntla et al. [38] as well as Gray and Ancey [17]. Extensions concerning density segregation as well as multicomponent segregation were already presented in the work of Marks et al. [29]. Other models for gravity-driven segregation mechanisms were presented by Larcher and Jenkins [25,26] or Hill and Fan [20]. An extensive summary about segregation in experiment and the mathematical models is given by Gray [15].
In this work, the focus lies on dry granular shear flows of systems of small and large particles, like they often appear in industrial and geological applications. In the upcoming section, a modified version of the hydrodynamic granular flow model of Latz and Schmidt [27] is introduced. The presented model additionally includes changes that consider the difference in particle size. Afterwards, a segregation equation is derived by means of mixture theory, similar to the aforementioned models. In contrast to the already mentioned models, a general direction of segregation is introduced. This aspect is relevant in order to make the segregation process invariant to the choice of the coordinate system. For this direction an expression is given, from which the segregation direction can be calculated for shear flows on the basis of the local flow field, and it does not have to be specified globally for the entire domain. In particular, this is done in three-dimensional space. Additionally, the segregation equation considers the compressibility of the granular flow model. The flux function of the segregation equation depends explicitly on variables of the granular flow equations. In Sect. 3, a numerical procedure is presented which considers the mentioned dependences of the flux function and prevents the system from reaching unphysical states. Further, the derived segregation equation reduces to the original model of Gray and Thornton [19] in case of parallel shear flows of an incompressible granular material. In the last Section, a first quantitative example is presented to show the interplay of the segregation equation and the modified granular flow model.

Mathematical model
The physical system of interest is given by a granular material of particles of different sizes and a passive air phase surrounding the material and filling the void spaces between the particles. Using the terms of the theory of mixtures, the whole system S consists of three constituents, the small particles S s , the large particles S l , and the air phase S a . It is assumed that the small and large particles of the granular material have equal material densities s * = l * = const . The volume fraction occupied by a constituent ∈ {s, l, a} is denoted by and it holds Consequently, the volume occupied by granular material is given by (2) s + l + a = 1.
In literature, it is often assumed that, e.g. in a granular avalanche, the volume of air a is constant. With the additional assumption of a passive air phase, it is argued that the air phase can be incorporated into the particle phases [16,18,19,37,39]. In such a case the mixture can be approximated by the particle phases which yields Nevertheless, in this work, it is assumed that the volume fraction of the air phase can vary locally. This is consistent with the compressible granular flow model which will be presented. For fulfilling (2), the air volume a must change according to the granular volume fraction . Certainly, the volume fraction of a particle phase relative to the granular material can be defined analogously to by Classically, equations are formulated in terms of partial densities . These are linked to intrinsic densities * , defining the quantity per unit volume of the pure phase. The linear volume fraction scaling is given by which equally holds for the granular mixture To avoid further variable changes, all equations are directly given in terms of the respective volume fractions, as all upcoming terms are linear in the density and * = s * = l * = const. The final model is given by a set of Navier-Stokes like equations for the granular mixture and additionally, an equation for the behaviour of the small particles as the large particle phase is given through the saturation condition (3).

The granular flow model
The model of choice to describe the behaviour of granular materials is a version of the Latz-Schmidt model [27], which originally describes a mono-disperse system of spherical particles. There are mainly two regimes that need to be taken into account when dealing with granular materials. First, the fast dilute flow regime which is dominated by binary particle collisions. The second regime is often called the static regime, where kinetic theory is no longer applicable, since the assumption of binary particle collisions breaks down for arresting granular flows. Therefore, the model is a hybrid combining properties of kinetic theory of granular gases, which is well verified in the field of dilute granular systems [4] and the framework of critical state plasticity which is employed in the field of (4) s + l = 1.
soil mechanics [36]. The model of Latz and Schmidt [27] has been validated for the case of inelastic, hard spheres against various experimental publications. A further extension to the kinetic theories has been given by Zemerli [42] presenting the transition of compacted, very dense flow via continuum mechanics of solids. This subsection serves as overview of the mentioned model as it has been derived and analysed in the already mentioned literature. Additionally, the model includes changes that consider the difference in particle size. The relevant expressions now additionally depend on the relative volume fraction ̂s . The general framework of the hydrodynamic model consists of three equations. The first Eq. (8) and the second Eq. (9) are the isothermal compressible viscous Navier-Stokes equations. They are solved for the granular volume fraction and the momentum , where defines the bulk velocity of the granular material.
The force term in the momentum balance (9) is assumed to be solely the gravitational force = .
Since it is assumed that the air phase plays a passive role, the interaction term is negligible, i.e. = . Further, p defines the pressure and ̄ the deviatoric stress tensor. In this model an asymmetric stress-strain relation is used, i.e., where is a viscosity. The usage of the asymmetric stress-strain relation is a simplification compared to the standard NSE. Typically, the symmetric version of the stress is used. In an d-dimensional space, it is given by In Eq. (12), is the bulk viscosity, and the identity matrix. The symmetric strain rate tensor ̄s is given by The choice of the less complex asymmetric stress-strain relation is not arbitrary. Asymmetric stresses are already used in the works of Mitarai et al. [31] and Campbell [5] to model collisional granular flow. Latz and Schmidt explained in [27] their decision to choose the non-symmetric tensor. For instance, they observed in numerical tests that the mixed derivatives in the symmetric tensor cause large unphysical spreading of granular jets perpendicular to the flow direction. Finally, they stated that at least in dense granular (10) = ̄, flow their approximation does not influence the results too strongly.
The third and last of the granular flow equations is the granular temperature Eq. (14). The concept of a granular temperature considers the energy transport due to particle movements and collisions. The granular temperature T, which is defined by the spatial average of the fluctuating part of the velocity, is a measure of the random particle motion in the granular system. One obeys the balance law where denotes the temperature conductivity. Equation (14) is quite similar to the classical heat transfer equation. It has a source in the viscous heating term 3 2 ̄∶̄ and a sink in the last term, where is the energy dissipation rate due to inelastic particle collisions. Please note, following a rigorous derivation as in [11] the heat flux includes not only a temperature but also a gradient in the number density and the dissipation term a gradient in the velocity. Here, the simplified approximation of the temperature equation given by [27] is used.
To close the system a constitutive relation for the pressure needs to be provided. As previously mentioned, one has to distinguish two regimes where the model should be applicable. Therefore, it is assumed that the closure for the pressure can be split in a purely kinetic and a so-called yield part These parts are subscribed by k and y, respectively. For the kinetic regime, Garzó and Dufty [11] derived expressions for the kinetic pressure and for the transport coefficients in the kinetic regime k , k , and k . Bocquet et al. [2] showed that there are simpler expressions than those derived by Garzó and Dufty [11] which produce quantitatively correct results in shearing experiments. Furthermore, Schmidt [35] showed that the simpler form together with the upcoming extensions is applicable to various regimes of granular flow. The simplified kinetic expressions for the transport coefficients are given by The compressibility factor (14) helps to ensure that the granular system stays in a physically valid state, i.e. < rcp (̂s) < 1 . The expression rcp (̂s) gives the upper bound for the volume of the granular system, named random close packing value. In contrast to the original mono-disperse model presented in [27], the random close packing value is not a constant any more, but it is a function depending on the local amount of small and large particles. An appropriate expression for the random close packing needs to be found by fitting experimental or simulation data. With the material parameters 0 , 0 , and 0 the granular material is characterized. It needs to be pointed out that the transport coefficients generally are a combination of species dependent transport coefficients. To keep the model as simple as possible, this first extension depends only through the radial distribution function on the species volume fraction to achieve nonconstant maximum packing values. Otherwise the amount of parameters which needs to be fitted increases or further relations need to be added to the model.
In the dense flow regime, the kinetic theory assumption of binary particle collisions is not valid any more. Due to the appearing long-term contacts, the presented relations (16) and (17) need to be modified. Similar as stated in the work of Savage [34], the pressure in the yield regime can be defined by The Heaviside step function ensures that the yield terms are only active in the dense flow regime, where the volume fraction of the material is larger than the cross-over volume, i.e. > co (̂s) . The cross-over volume is smaller than the random close packing volume, but it is assumed that it behaves the same way depending on the local amount of small and large particles. It is close to the random loose packing value. The yield parameter T 0 is a positive constant, it prevents the pressure to vanish for zero granular temperature. Similar to the pressure, the transport coefficients have to be modified in the dense regime. The final terms are defined in the work of Schmidt [35] Due to definition (20) the transport coefficients increase for diminishing granular temperature. This behaviour can be observed in experiments [9]. Although the form of the constitutive relations is held relatively simple, the model reproduces properties of granular systems like Bagnold scaling, frictional properties, and plastic regime as it has been shown by Zemerli [42]. The granular material is uniquely characterized by the macroscopic parameters rcp , co , 0 , 0 , 0 and T 0 . These parameters generally have to be validated by experiments, but in many cases, some parameters can also be approximated via further relations. For instance, Latz and Schmidt [27] used the relation to determine one of the coefficients 0 and 0 with the help of the internal friction coefficient .

Segregation modelling
To derive a segregation equation that fits the granular flow model, the modelling process of Gray and Thornton [19] suits as orientation. Based on the assumptions of mixture theory, it is assumed that for each phase a mass and a momentum balance are given by The exact form of the several terms in (22) depends on the assumptions made for the physical system.
First, it is assumed that the particles in the granular system do not amalgamate or break. This implies that the mass transfer variable is equal to zero, m = 0 ∀ ∈ [s, l, a] . Consistently, as stated for the granular flow model, the body force on each phase is solely the gravitational force, = ∀ ∈ [s, l, a] . The stress tensors for the particle phases in the segregation process are pressure dominated. They are approximated by lithostatic pressure fields, i.e. ̄ = −p [19]. The assumption of a passive air phase implies that the amount of stress in the mixture coming from air is negligibly small, i.e. ̄a = , and also the interaction force a = . Hence, one can focus on the particle phases and the combined granular system as it is usually done in literature.
With the assumptions made so far, the balance laws (22) for a particle phase are given by (21) The assumptions imply that the pressure of the whole mixture is equal to the granular pressure, i.e.
As previously mentioned, there is no interaction force between the particle phases and the air phase. Hence, An equation describing the change of one particle phase can be derived, starting from balance law (23a). Introducing the relative velocity between phase and the bulk the mass balance can be rewritten in the form In this form, the left-hand side of Eq. (27) describes the transport of particle phase with the bulk velocity . The right-hand side represents the segregation. It is a motion relative to the bulk with velocity ̃ . Segregation is a result of the interaction between the small and the large particles. This interaction is given in the momentum balance due to the interaction force . Hence, an expression for the relative velocity can be obtained from the momentum balance (23b).
With the assumption that the partial densities and momenta become quasi-steady even before the segregation starts, like in [39], the momentum balance (23b) reduces to Summing (28) over all constituents implies that the pressure field is lithostatic, The key idea to derive an expression for the relative velocity is to introduce a pressure scaling that differs from the standard volume fraction scaling as it is used for the densities. The scaling, where the partial pressure depends linearly on the bulk pressure was formulated by Gray and Thornton [19]. The idea arises from the assumption that small particles carry less of the overburden pressure than the large particles while they percolate to the ground of the mixture. The factor f determines the proportion of the load carried by phase . Equation (30) automatically implies that If solely one of the particle phases is present, it has to carry all of the granular load. Hence, the proportionality factor f must fulfil condition Similar to [19] the expressions for each proportionality factor are chosen in such a way that they satisfy the conditions (31) and (32), i.e.
The expressions in (33) are not unique. In this work, this simple form is used, where B is a non-dimensional magnitude. It should be mentioned, that Gajjar and Gray [10] discussed different expressions in their work, where they modelled asymmetric flux functions. Furthermore, an expression for the interaction force is required. This expression is modelled by which satisfies the summation condition (25). The last term of expression (34) is a Darcy term, where c is the coefficient of inter-particle drag. Its nature is poorly understood, thus it has been chosen to be constant in the first works modelling segregation in that way [18,19]. The Darcy term was already introduced by Morland [32] after observations that the segregation process of particles shows an analogy to the percolation of fluids through porous solids. The first term ensures that the percolation is driven by intrinsic than partial pressure gradients as in Darcy's law [19]. The segregation depends on the shear rate ̇ . It is inversely proportional to the drag term as an increase in the shear rate decreases the friction between the particle layers. The first time it appeared in a gravity-driven segregation model, was in the work of Marks et al. [29]. It is allowed to vary in this work depending on the flow field. The expression for the interaction term (34) is the established form in literature when dealing with gravity-driven segregation models [10,16,19,38].
With the definition of the interaction term (34) and the pressure scaling (30), Eq. (28) is given by In the most works (e.g. [3,6,8,16,18,19,25]) a segregation equation is derived for avalanches flowing down inclined planes with a constant angle. Further, it is assumed that the avalanche is incompressible with a constant height. In these gravity-driven shear flows, shear bands appear that are parallel to the inclined plane. In such a case, the coordinate system can be rotated such that the avalanche flows in the positive x-direction and the segregation solely takes place in the z-direction. Especially, this means that the segregation direction will not change, neither over time nor locally in the spatial domain. To model the segregation process in such a case, one can focus solely on the z-direction as it has been done in the already mentioned literature. In this work, the segregation equation should be coupled to the granular flow model where arbitrary flow directions are possible. Therefore, the modelling done here should be more general. Assume that the segregation direction at a point in space is given by the vector s . This vector spans a onedimensional subspace of ℝ d , denoted by D . To focus on this subspace one can define the projection ∶ ℝ d → D by To be more precise, the direction of the segregation is perpendicular to the shear layer which depends on the velocity field.
Applying the projection on Eq. (35) gives due to the linearity of . As the pressure field is lithostatic (29), Eq. (37) rearranges to Note that the relative velocity ̃ is the segregation velocity. It lies in the image of the projection function and therefore per definition it holds Finally, using the expressions for the proportionality factors (33) the relative velocities for the small and large particle phases are given by With the expressions for the relative velocities (40), the balance law (27) modifies to respectively. The expression has the dimension of a velocity and gives a measure for the segregation speed of small particles in the system. It holds ls = − sl . As previously mentioned, the shear rate is assumed to be non-constant in this work. It should depend on the state of the granular flow model. This can be done using the granular temperature. The shear rate can be defined by With definition (43), the viscous heating term in the granular temperature equation (14) is simply given by 3 2̇2 . The shear rate in the granular system can be approximated by the granular temperature. In a stationary uniform state, the equation for the granular temperature (14) can be written as Here, it is assumed that the heat flux term is small compared to the viscous heating and dissipation. Then, the temperature equation is in the equilibrium between viscous heating and dissipation [42]. With the expressions given in (17) and (20), the shear rate is proportional to the square root of the granular temperature With the derived expression (45), the segregation velocity can be rewritten as where the segregation rate S sl 0 is constant. In general, the segregation rate depends on the properties of the respective granular material. Therefore, the segregation rate needs to be fitted to the real granular system. Furthermore, S sl 0 is expected to be positive. Finally, the segregation equation for the small particle phase takes the form The structure of the derived segregation equation (47) is similar to the models mentioned in Sect. 1. The main differences are the non-constant granular volume fraction and the dependence on the granular temperature T, which is non-constant as well. Not least, this will change the way the equation has to be handled numerically. For simplicity, the superscript of the volume fraction of the small particle phase is dropped = s , since only this phase is used in the following. One should note, assuming constant granular temperature and constant granular volume fraction in the granular system the model reduces to the original model of Gray and Thornton [19]. For the sake of completeness, it should be pointed out that the dependence of segregation on granular temperature has already been discussed in the work of Fan and Hill [8]. The dependence arises from the assumption that segregation is driven by kinetic stress in the definition of the interaction term of the momentum balance.

The segregation direction
For the modelling of the segregation equation, it has been assumed that a vector s is given that spans a subspace of ℝ d in which the segregation takes place. To solve the segregation equation such a vector needs to be found. It is well known that for gravity-driven shear flows, the segregation happens in the direction perpendicular to the local shear layer [3,19,26]. The shear depends on the velocity field of the granular system or more precisely on the gradient of the velocity ∇ , such that s = s (∇ ).
To find the segregation direction in granular shear flows, a concept from the field of mechanics is adapted. Assume a mechanical body under some stress. This stress state is given by a symmetric stress tensor ̄ with principle stresses and corresponding normalized eigenvectors 1 , 2 , and 3 , which are mutually orthogonal. As described, for example, in the work of Wu [41], planes of the principle shear stresses ij of the body are orientated as depicted in Fig. 1. The plane of maximal shear stress is defined by the eigenvectors corresponding to the smallest and largest eigenvalue. Under assumption (48) the maximal shear stress max = 13 acts in the sectional plane depicted in Fig. 1c.
Adapting this procedure to the framework of granular flows, the segregation direction can be computed by using the symmetric versions of the stress and the strain as given in (12) and (13). From the stress-strain relation, it can be deduced that the eigenvectors of stress and strain are equal. The corresponding eigenvalues might be different but have the same order in size. Hence, the planes of principal shear stress and principal shear strain coincide. As solely the orientation of the shear planes is of interest, it is enough to look at the symmetrized velocity gradient of the bulk flow given in (13) by Unfortunately, normalized eigenvectors are not unique with respect to the algebraic sign. Hence, the planes depicted in Fig. 1 are not unique.
For mechanical bodies, the forces acting in the defined planes between, for example, all combinations of ± 1 and ± 2 are equal (see Fig. 2). This is not the case for the granular system as it is not one solid body. Hence, it is assumed that the shear plane is the one which is mostly parallel to the velocity field. Hence, the following procedure can be stated to compute the segregation direction in granular shear flows: 1. Compute the eigenvalue max and min of ̄s = 1 2 ∇ + (∇ ) T . 2. Compute the corresponding eigenvectors max and min . 3. Define s = max ± min depending which is mostly orthogonal to .

Numerical treatment
The computation procedure for the previously presented set of equations can be done in the following way. From some initial state of the granular system the local values for the random close packing rcp (̂) and the cross-over volume co (̂) can be computed. Then, the granular flow equations are solved to gain the granular velocity and the scalar fields for granular volume fraction and granular temperature T. With these data the segregation equation is solved to update the distribution of the small particles in the granular system. This process can be repeated for each time step starting again with updating the new local values of rcp (̂) and co (̂) from the computed -field. The granular flow model can be solved using a finite volume approach as described in [35]. Due to the compressibility of the granular system one needs to be more cautious solving the segregation equation. Whereas the transport with the granular bulk material should be solved similar to the mass balance of the granular flow system, the explicit dependence of the flux function on the granular volume fraction needs to be treated in a special way. Since ∈ [0, ] and = (z) is varying locally, a suitable numerical scheme needs to be chosen to prevent the system from reaching unphysical states, where overshoots . One has to tackle the problem of a spatially dependent flux function.

Spatially dependent flux functions
Since the segregation flux is the point of interest deriving a suitable numerical scheme, the segregation equation (47) is reduced to the relevant parts given by a one-dimensional equation of the form For now, it is assumed that the z-direction is the direction of segregation. Hence, the bulk flow can be neglected as it is assumed to be perpendicular to the segregation. For simplicity, the granular temperature is assumed to be constant. The spatial dependence is given only due to the granular volume fraction.
As stated by LeVeque [28], the spatially dependent flux function can be discretized and one obtains a flux function = 0. f i ( ) associated with the i-th grid cell which yields a Riemann problem at z i−1∕2 with states n i−1 and n i by Solving such a Riemann problem generally consists of finding a state ̌L that can be connected with left going waves to n i−1 and a state ̌R that can be connected with right going waves to n i , such that additionally One or both of the states ̌L and ̌R do not have to coincide with the states n i−1 and n i . A flux approximation obtained from a classical scheme (e.g. Upwind) of the directly discretized flux fails in this case, because the chosen states do not fulfil condition (51).

A computational method from Riemann solutions
To find a numerical scheme fulfilling (51) can be done by finding the solutions of the appearing Riemann problems. As presented by Jin and Zhang [23] solving a problem in the field of traffic flow, one can extend Eq. (49) to a nonlinear resonant system of equations, as defined in [21], by introducing an additional balance law for the spatially dependent variable of the flux term.
Since bulk transport and segregation act orthogonally to each other, Eq. (52b) is a suitable approximation of the mass balance in the spatial z-direction. The linearised version of the upper system is explicitly given by The eigenvalues of the Jacobian are and the corresponding right eigenvectors are (51) (55) From the eigenvalues, it follows that the system is nonstrictly hyperbolic since, for several states in the system, it holds A state * = ( * , * ) T is called critical if In the -space, all critical points form a smooth curve. This curve is called transition curve and is defined by A Riemann problem is given by Eq. (52) and the initial condition Associated with each eigenvalue, one basic wave solution exists, travelling with speed 0 and 1 , respectively. The integral curves, the solutions travel along, can be constructed from the eigenvectors. For the given problem, they are given by all states where f ( ) = const for the standing wave (contact discontinuity) with speed 0 = 0 and all states where = const for a wave travelling with speed 1 . The respective integral curves passing through a critical state * and the transition line are depicted in Fig. 3.
Isaacson and Temple [21] found that resonant systems, as derived above, can be uniquely solved by introducing a new (57) 0 = 1 . Fig. 3 Plot of the integral curves passing through a critical point * = ( * , * ) T and the transition line in state space entropy condition, which states that standing waves must not cross the transition line (59) in state space. Similar to [23] and depending on the location of the left Riemann state L in relation to the transition line, one can find ten different combinations of waves to connect L to a state R . Hence, there are ten different regions A-J in the state space for the position of R that lead to one of the ten wave combinations. In Fig. 4, these regions are plotted. Six regions if L is located left of the transition line and four if L is located right of .
Since the volume fraction of a particle phase must not be larger than the granular volume fraction , all states in the state space where > are forbidden. All solutions are a combination of two or three different waves, where each wave is either a shock wave S , a rarefaction wave R , or a contact discontinuity C . An example of a solution, consisting of three waves, where L is located left of the transition line and R in region is shown in Fig. 5. The correct flux approximation, fulfilling (51), for each of the ten different cases is given in Table φ-axis  Table 1 and the flux function (63) are comparable with those given in [23] but for a convex not a concave flux function. Therefore, the notation of [23] has been adapted.
As derived in Sect. 2.2, the spatial dependence of the segregation equation is not only given due to the granular volume fraction but also due to the granular temperature T. Solving the equation numerically, this additional dependence makes no big difference. Analogously to the statements previously done, the numerical scheme can be used almost in the same way including the granular temperature [14]. Hence, the explicit description of the flux function can be written in the exact same way as in the constant temperature case, where this time Finally, the multidimensional version of the segregation equation can be solved as series of one-dimensional equations using Godunov type dimensional splitting. Constructing a regular mesh of finite volumes, the discretized

Computational examples
The previously derived method is used to simulate the mixing of a granular material in a rotating tumbler, but first of all one simple one-dimensional configuration serves as test case to show the advantage of the previously derived numerical scheme.

Gaining physical solutions by the modified Godunov method
Assume the simplified version of the segregation equation (49) with non-constant granular volume fraction and also with non-constant granular temperature T given by The spatial domain used, is the unit interval [0, 1] in z-direction, which is assumed to be the direction of the segregation process. For the simulation, a spatial grid of 50 cells is used and the time step satisfies the Courant-Friedrichs-Lewy (CFL) condition. The segregation rate is chosen to be S = 1 . The simulation result after 2 s is shown in Fig. 6. The granular volume fraction and the distribution of the granular temperature are chosen arbitrary, depicted by a dashed and a dotted line, respectively. Initially, the system is in a perfectly mixed state where = 2 . One can see that the classic Godunov or Upwind scheme, where simply the exact value for and T at the cell interfaces is used for the flux evaluation, produces unphysical states, because relation (51) is not fulfilled. In [14], it has been shown that the overshooting error for the classical schemes decreases for decreasing grid size, but never vanishes. In contrary, the modified Godunov scheme always holds the system in a physical state, where ≤ .
In a close up view the relevant part at the presented time is depicted again.

The three-dimensional rotating tumbler
The chosen case to test the performance of the whole model is a thin, plain rotating tumbler with diameter 0.1 m. It is discretized with 100 finite volumes per diameter and periodic boundary conditions at the end walls. This guarantees a large granular bed and lower computational costs compared to simulating a tumbler of extended thickness. For the granular material, glass particles are used similar as in [33]. The used parameters to characterize the granular material are given in Table 2.
For the functional parameters of random close packing rcp and the cross-over volume co a polynomial fit is used, based on the data presented in [24] and the segregation rate is set to S sl 0 = 1. Initially, the system is perfectly mixed and the tumbler is filled to 75 % as shown in Fig. 7a. To simulate the rotation of the tumbler the whole set of equations is solved in the rotating reference frame. For the illustration, the whole grid is simply transformed back to the inertial system. The rotation speed is set to 0.1 revolutions per second (rps). The segregation direction for each cell is computed directly from the velocity field, as presented in Sect. 2.3. Figure 7b shows the result for the computed directions in a close up view of the surface region of the granular system. Coloured arrows depict the velocity field and black ones the computed segregation direction. As expected for avalanche-like shear flows, the segregation directions are orthogonal to the velocity field.
The general behaviour of the granular material in such a tumbler is already well-known, since it is used in industrial processes and therefore, a thoroughly investigated case in literature. Furthermore, the granular material shows different regimes during the rotation of the tumbler. There are flowing regimes where segregation takes place as well as static regimes.  Fig. 6 a Simulation result of the one-dimensional segregation process in z-direction for the small particle phase, given by Eq. (66) at time t = 0.75. The black dashed line depicts the granular volume fraction , which is an upper bound for the small particle phase. The distribution of the granular temperature is shown by a dotted line and the segregation rate is chosen to be S = 1. Initially, the system is perfectly mixed, which is depicted by a light grey line. Over time, the small particles concentrate at the left. Only the modified Godunov scheme prevents the particle phase from overshooting . b A close-up view shows the problematic region  The behaviour of the granular material in a rotating tumbler depends not only on the material properties but also on the rotation speed and the filling height. Mellmann [30] gives an extensive overview of the different conditions and regimes.
For the chosen test case a ring-like formation of small particles, large particles near the outer walls, and a mainly untouched core of the granular bed are expected. Figure 7c depicts the simulation result after 25 s of rotation, showing the expected segregation pattern.
Due to the modification of the random close packing rcp (̂s) and the cross over volume co (̂s) compared to the original model in [27], the formation of the segregation pattern also influences the volume fraction profile of the granular bed. For a mono-disperse particle system, as assumed in the original version, the volume fraction increases from the surface to the bottom of the granular system. Hence, the volume fraction distribution is similar to a lithostatic pressure profile. In case of a binary particle system, now this lithostatic profile is superimposed by a profile that is similar to the segregation pattern. One can see in Fig. 7d that the volume fraction is much higher in regions where the system is partially mixed than in regions where solely one particle phase is present.
The presented results in Fig. 7 can only serve as qualitative test to show the general behaviour of the combined model. For a detailed comparison with experimental data, more effort has to be put in the investigation of the characterization of the granular model. As shown in [14], the computed fields of the granular flow model have a strong influence on the segregation equation, specifically the distribution of the granular temperature. Therefore, it is important to find a suitable parameter set for the granular flow model, that accurately characterizes the real granular material.

Conclusion
This paper investigates the simulation of granular material flows with segregation. The presented model extends a set of Navier-Stokes-like equations with a non-Newtonian and non-linear rheology model, to incorporate the segregation process of binary granular systems with small and large particles in up to three space dimensions. The segregation equation has been derived for shear flow applications using mixture theory. The segregation process is regulated due to the dependence on the granular temperature. The model, as derived in this paper, depends explicitly on the local value of the granular volume fraction. The volume fraction is non-constant over the whole domain, since the granular flow equations are compressible. This entails that special care needs to be taken in its numerical treatment. Therefore, a modified version of Godunov's method has been formulated for solving the segregation equation, preventing the system from leaving a physical state. To make the model applicable to a wide field of applications, where the flow direction might change, a local segregation direction vector has been introduced in the three-dimensional space. This segregation direction is computed based solely on the granular flow field, which also makes the segregation process invariant to the choice of the coordinate system. An advantage is that the segregation process does not have to act in a single direction in the whole domain. First simulations with the combined model were done with the test case of a rotating tumbler. The segregation pattern observed and the volume fraction profile for the case of a binary particle system are consistent and fulfil the expectations.
In this elaboration, the model is held as simple as possible for the first qualitative tests. The starting point is an solely advective equation, which is limited to the topic of size segregation. The authors are aware of the fact that more progress has already been made concerning diffusion terms, more complex flux functions, or the dependence on shear rate and granular temperature. However, extending the model by a diffusive term or a combination of size and density segregation, as already done in literature (e.g. [17,18,38]), should be possible. Instead, this work focused on the topics compressibility, the segregation direction, and the numerics for physical solutions in the field of granular shear flows.
Note that for applications in a kinetic regime, several terms, like the transport coefficients in the granular flow equations need further modifications. These can be found by Enskog theory for polydisperse mixtures, as presented by Garzó et al. [12,13]. The same holds true for the segregation equation. Leaving the dense shear flow regime, modifications in the segregation terms seem to be necessary.
Furthermore, in the dense regime, the set of granular flow equations already has a handful of parameters which characterize the granular material. Hence, a detailed analysis of the model and comparison with experiments in further case studies will be relevant for a better understanding of the interplay of segregation and granular flow equations. This is planed in a connecting elaboration.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.