Flows of granular material in two-dimensional channels

Secondary cone-type crushing machines are an important part of the aggregate production process. These devices process roughly crushed material into aggregate of greater consistency and homogeneity. We apply a continuum model for granular materials (Jop et al., Nature 441:727–730, 2006) to flows of granular material in representative two-dimensional channels, applying a cyclic applied crushing stress in lieu of a moving boundary. Using finite element methods, we solve a sequence of quasi-steady fluid problems within the framework of a pressure-dependent particle size problem in time. Upon approximating output quantity and particle size, we adjust the frequency and strength of the crushing stroke to assess their impacts on the output.


Introduction
In this article, we will model flows of granular material in cone-type rock crushers, as depicted in Fig. 1. A crusher chamber typically consists of an inner 'mantle' that is driven by an eccentric shaft within an outer static 'concave', as idealised in Fig. 2. We will utilise the viscoplastic continuum relation for granular materials devised by Jop et al. [1], and approximate two-dimensional flows of pressure-crushable granular media within channel cross sections similar to that shown in Fig. 3. Material to be crushed falls between the two surfaces, compacting locally where the inner mantle has the greatest proximity to the concave. The mantle is fixed at its summit, while its base rotates on a circular path about the concave's axis of symmetry, as in Fig. 4. We will consider the channel to have length L and height h at the inlet. The outlet height h out at a given point on the circumference of the crusher will oscillate between the minimum and the maximum sizes at between four and ten cycles per second. Within the literature, the most prominent existing model for material flow and crushing within crushers was developed by Evertsson in his doctoral work [2,3] from the earlier population model [4]. That approach was to establish a mass balance in combination with a flow and crushing model. The flow model is based on computing the path of each single particle through the crushing chamber, whether they be in freefall, sliding or crushing. The system is controlled by a 'volumetric filling ratio', η v < 1, which allows for the incomplete filling in the inlet zone seen in real crushers. This model was then followed up by Evertsson and Lindqvist in collaboration with others [5][6][7][8][9][10][11] with the main intentions of improving the flow model and integrating lining wear into simulations. They optimised the geometry such that material experienced no sliding at the basal surface in its transit through the cross section. In an independent succeeding study [12][13][14][15], Evertsson's results were refined further and successfully applied to an operational crusher and are now being incorporated into new crusher design. The Evertsson model is interesting due to its discrete approach to the system, but the ultimate integral equation is a far less tractable model than our continuum fluid flow approximation, as we aim to examine the overall behaviour of the material and the quality of the output rather than to simulate the explicit grain trajectories and crushing.
The modelling of granular media has become increasingly popular in recent years as researchers attempt to understand the three main modes of behaviour-which equate roughly to gaseous, liquid and solid regimesand generate the tractable, unifying models that still seem to be beyond our grasp. Much time has been spent into modelling the gaseous phase using 'granular temperatures' [16], and the solid phase via solid mechanics [17], but the dense liquid flow has defied easy modelling. The dynamics/mechanics involved are complex, given the irregularly shaped particles which arise within any natural granular material, and random assemblies of these particles can lock together into temporary structures, which we refer to as force chains or networks [18][19][20][21]; which structures are sometimes also referred to as force eddies or clusters [22].
Instead of a discrete model, we choose to adopt a contiuum model that is minimal in nature and is based empirically on the study of the pseudo-liquid regime of flows of dense granular media, on an assumption of rock particles being small enough for a continuum approximation. The complexity and lack of consensus on modelling granular materials in the literature also make such an approximation credible in the context of our problem. Our chosen model is in essence composed of a pressure-dependent yield stress and an effective viscosity depending on the pressure p, grain diameter λ and shear rate |γ |. In common with what we assume to occur inside the crushing cavity, this material is pressure-thickening, which is analogous to the slowing and cessation of the fragments in the crusher as the crushing stroke occurs. Employing techniques used in computations involving Bingham viscoplastics [23], we choose to regularise the dimensional constitutive using the Bercovier-Engelman method. We can then apply the formerly piecewise and now universal constitutive relation over the whole channel Ω. The choice of regularisation is arbitrary in the sense that we could choose the Papanastasiou [24,25] method equally well in the context of this work. With the prerequisite analysis, we might invoke variational inequalities [26,27] in place of regularisation, prompting a radically different formulation. We choose to invoke grain size dependence, in this instance purely by its directly stated relation to the inertial number I : where deviatoric stress τ is given by and total stress σ by for close-packed material density ρ s , the rate of strain tensorγ , effective viscosity function μ(I ), and identity tensor I . In reality, many of the parameters, such as the coefficients ζ s , ζ 2 , inertial constant I 0 , packing fraction κ and Young's Modulus E c might also depend on λ. For the purposes of this computation and simplicity, we restrain dependence as stated, and neglect the influence of polydispersity in grain shape. Further λ-dependence could be investigated in future work.
In the main body of this paper, we will apply slip boundary conditions at both the physical boundaries. In the preliminary unidirectional work in Sect. 2, we prescribe no-slip at the basal surface for a unique solution in the case of a parallel-sided channel. We prescribe basal slip in the ensuing two-dimensional work due to the assertions in [3], wherein it is stated clearly that there are three modes of movement for material within cone crushing machines: free fall from the upper surface, ascendancy on the basal surface with the crushing concave and sliding on the basal  [30] surface. Evertsson states slip is common [3] but then optimises his model to avoid slippage. Within [28], it is again stated that slip is common within the first few layers of grains at the base. We proceed by presenting a preliminary set of unidirectional solutions in Sect. 2 which demonstrate our slip boundary condition and regularisation before moving on to the two-dimensional work. The variational formulation for the problem is given in Sect. 3, the procedure for performing the calculations in Sect. 4, some illustrative results are given in Sect. 5 and a discussion on the two-dimensional system is contained in Sect. 6.

Preliminaries
For ease of reference, we list typical dimensional values for quantities involved in Table 1.
To begin, we examine the steady, non-inertial, gravity-driven and incompressible flow of granular material in parallel-sided channels inclined at an angle θ to the horizontal, as shown in Fig. 5.
From our Table 2 of typical nondimensional parameter values, we can deduce an approximate Froude number F n of the form which number can stand in as an estimate for the Reynolds number in this case due to the pressure dependence of the viscous stresses. Since F n ≈ 10 −1 , we can consider the flow to be broadly non-inertial.  Unidirectionality dictates that and p simplify to functions only of y, the distance from the channel base. We nondimensionalise the system using the scalings √ gh for velocity u, ρgh for stress-like quantities τ , p, E e , h for coordinates x = (x, y) and grain diameter λ, √ h/g for times t and T c , and gh 5 for volumetric flux Q. ρ denotes material density, and g denotes gravitational acceleration. In line with this, the effective viscosity μ [1] simplifies to the following regularised and unregularised forms where is a small non-negative regularising parameter applied in the Bercovier-Engelman fashion [23], and the deviatoric stress is then defined by Our aim in preliminary studies is to find solutions of the unregularised and regularised JFP materials which may then be used to illustrate and validate our choices of slip condition and regularisation. These solutions may also be used as initial conditions for the two-dimensional calculations if necessary.
As in most fluid problems, we will solve the system of equations formed by the equations of conservation of mass and momentum. Under our unidirectional assumption, conservation of mass is described by and is automatically satisfied, while conservation of momentum simplifies to is the vector of gravity forces. We prescribe the boundary conditions: and employ the specific slip relation from [31], wherein tangential deviatoric stress τ nt and velocity u t vectors are defined at the wall by and n and t denote normal and tangential vectors to the boundary, respectively. We prescribe the boundary stresses as both velocity and pressure dependent via The quantity ζ w is the material friction coefficient. ζ s , ζ w , and ζ 2 can be alternatively understood via their relations to friction angles φ s , φ w , φ 2 , The quantities τ ω and u ω are our boundary condition smoothing constants and represent small non-negative stresses and speeds, respectively, such that the approximation collapses down to the piecewise slip condition:  Convergence plots for τ h for various configurations of (u ω , τ ω ). In this case, θ = 55 We in effect prescribe that for a material to slip at the boundary, it must reach a threshold shear stress τ slip ∝ p, thus making slip less likely for higher applied crushing pressure p h . The term (τ ω /u ω )u t exists as a guarantee of a bijective relation between slip stress and slip speed. Within the preliminary work, it is simply implemented, while in the two-dimensional finite element calculations, it is prescribed as a natural boundary condition, which requires the addition of complementary zero-flux terms in the variational formulation. For other examples of slip boundary conditions and discussion see [28,31,32], while for a variational inequality possibility as applied to the Signorini boundary condition problem see [33]. In this instance, we choose regularisation in line with our previous choice in the effective viscosity function μ(I ). Under unidirectionality, our prescribed slip stress τ slip (4), simplifies to the scalar: which, in the piecewise case, would dictate where u h = u(h) is the slip speed. As a verification of our smoothed boundary condition, we have run a sequence of calculations, and display the results in Figs. 6 and 7. We do indeed see that the smoothed boundary condition converges to the piecewise condition with the decreasing u ω and τ ω /u ω in both the speed and stress fields. We will solve the unregularised JFP system analytically with the exception of one key parameter y 0 which must be found numerically in the event of nonslippage, while the regularised JFP system will be solved numerically using finite difference techniques.
where s = sign (du/dy). Since the material is unregularised, we expect the velocity profile will have a rigid plug bounded by y 1 , y 2 and consisting of at most three distinct regions as shown in Fig. 8. The constant y 0 determines the zero point in the linear stress function, and is used to indirectly prescribe the boundary condition on the stress at the surface y = 1. The number of regions present is not constant and depends primarily on the applied normal stress p h . For the special case of a free surface situation when p h = 0, we anticipate that the material flow will consist of a single shear layer, but for all p h > 0, an unyielded region should exist. In general, we expect the following three modes of behaviour to occur as we vary p h : -Free surface flow for zero p h , one shear layer, -Slip for intermediate p h , a plugged region and one or two shear layers, -No-slip for large p h , a plugged region and two shear layers.
We solve Eq. (7) with the boundary conditions To integrate, we use the substitution: where, under this substitution, The equation for u(y) may be simplified as follows: We then integrate for the general solution: and apply our boundary conditions to find a piecewise solution. For 0 ≤ y ≤ y 1 , we have while for y 1 < y < y 2 , we find u(y) = u(y 1 ), and finally for y y ≤ y ≤ 1, we have It should be noted that the values of A s , A 2 are dependent on s, and so have distinct values in the two possible shear layers.
To determine the plugged region boundaries y 1 , y 2 , as shown in Fig. 8, we solve du/dy = 0, from Eq. (7): The boundary conditions on u(y) ensure continuity of the solution at y 1 , y 2 . The shear stress at the surface τ h is determined indirectly by the chosen value of y 0 . To find a solution consistent with the stress boundary condition, we must numerically find a root in y 0 of the function: The surface velocity u h is determined by modifying the normal applied stress p h or the friction coefficient ζ w , but this may not be done arbitrarily as there are bounds on these two quantities which we shall now derive. Under unidirectionality the stress-strain relation for an unregularised JFP material, as detailed in [1], is given by , and therefore Since, by the formulation of the constitutive relation in [1], we must have it follows that at the base and then at the upper surface, We combine the two inequalities on the assumption that there is an overlap in the regions for y 0 , and hence and gain the final requirement for p h : The bound on p h can be physically understood as the limiting pressure for material flow. Above this limit, the material is locked into a static arrangement due to its interlocking grains. The left-hand side of (11) is uniformly negative for θ ∈ (φ s , φ 2 ), and so practically we have since non-negative (pulling) applied normal stresses are nonphysical in this case. Further conditions on parameter ranges can now be determined based on the physical situations we wish to model and practicalities. For free surface flows, we will need p h = 0 and therefore a negative lower bound In the more general case of flows where p h > 0, we look for other constraints; we require a positive upper bound on p h , Steady flow is therefore considered realistic for the JFP model under a bounded applied normal stress p h if while no flow is possible for tan(θ ) < ζ s for whatever p h , and there speculatively could be unsteady flow for tan(θ ) > ζ 2 for some values of p h . In light of these bounds on applied stress p h (12) and channel inclination θ (13), we modify the original values of ζ s , ζ 2 . For sufficiently large p h , we require small ζ s , and for large θ , we require large ζ 2 . Limits for the the slip coefficient ζ w are also built into the model. Since the edges of the unyielded 'plug' region are given by Eqs. 9a, 9b, we can deduce that for the existence of a second shear layer and hence 'controlled slippage', we require The ratio ζ w /ζ s is, in addition to p h , vital for the onset of slippage as ζ w /ζ s 1. We may also speculate a similar change in behaviour with the ratio ζ w /ζ 2 where slip decreases as ζ w /ζ 2 → 1 for given p h . Given all of the above, we can plot solutions such as those shown in Figs. 9 and 10. It is apparent that our slip condition is operating effectively as slippage increases with decreasing applied stress p h . The slip velocity u h approaches but does not cross zero with increasing p h . We interpret positive u h as unphysical in this case.

The regularised JFP system
We regularise the JFP constitutive relation and introduce a parameter > 0 such that the effective viscosity is given by (1). With the inclusion of , we no longer have a tractable analytical problem and resort to numerical solution via finite difference techniques. We discretise u, τ and p over y as in Fig. 11, and approximate the governing equations: and boundary conditions by the following system of 2n − 1 equations This is then solved using Newton's method. We demonstrate convergence of the regularised solution to the unregularised with decreasing in Figs. 12 and 13. The regularisation effect is to increasingly eliminate the plugged region in the velocity profile that we see in the unregularised profiles in Fig. 9. The regularisation is especially apparent in the stress profile where the discontinuity is clearly smeared, illustrating how over-regularisation fundamentally distorts the nature of the solutions computed. The regularisation parameter may be characterised as a function of the Jump number [34], a quantitative measure of the difference between the yielded and unyielded behaviours.

System
We now move on to examining flows of two-dimensional granular materials within nonparallel-sided channels, such as that shown in Fig. 14. The geometry function h(x) is given by which has been designed to satisfy conditions We retain our previous scalings and our non-inertial and incompressible flow assumptions. The whole timedependent problem is modelled as a sequence of quasi-steady fluid problems. Since this viscoplastic constitutive  (14), and basal service y = 0. Applied stress p h is prescribed at y = h(x) relation was designed for generic granular materials, we modify some of the model parameters for quantitative similarity to our real-life crusher. The crusher's cross section is approximated as a two-dimensional channel of fixed geometry. The fixed geometry is chosen since the moving boundary problem would be problematic in combination with our assumption of incompressibility, while the crushing stroke is simulated by a prescribed applied stress p h (t) at the inlet. This idealised cross section is a minimal model that represents approximately the effective crushing portion of the cross section. The pressure drop from inlet to outlet is a necessary compromise to drive the fluid continuum approximation to the granular material, in addition to the gravity forces which will cause flow under sufficiently small applied normal stress.
We will solve a sequence of quasistatic problems in velocity u(x; t) and pressure p(x; t) for fixed λ in space, linked by a constitutive problem in time and space for λ and fixed u, p. Assuming small enough time steps, we can assemble the solutions to these two problems into a larger composite solution for (u, p, λ).
The two dimensional viscoplastic fluid problem in velocity u = (u, v) and isotropic pressure p is composed of the equations of conservation of mass and momentum as follows: which are equivalent to the following in terms of u = (u, v) and p, where we have the regularised nondimensional constitutive relation defining the effective viscosity φ as nonlinearly dependent onγ , p and particle size λ by φ(|γ |, p) = ζ s ( 2 + |γ | 2 ) 1/2 + λκ 1/2 (ζ 2 − ζ s ) The volume fraction κ is assumed constant in line with our assumption of incompressibility. This system is accompanied by the following boundary conditions: u · n = 0 at y = h(x) and y = 0, (τ · n) · t = τ slip (u) at y = h(x) and y = 0, p = p in (y) at the inlet (x = L), p = p out (y) at the outlet (x = 0), v = 0 at the outlet and inlet, σ · n = p in (y)n at the inlet, σ · n = p out (y)n at the outlet, where the pressure p has been carefully prescribed at the inlet and outlet due to the pressure dependence in the constitutive relation. The pressure conditions at the inlet and outlet are completed by which describe oscillatory applied stresses at those boundaries with a finite pressure drop p step from inlet to outlet, bounds at the inlet for the applied stress of p h,low and p h,high , and crushing stroke frequency R.
The partial differential equation for the particle size λ is given by where λ = h out at t = 0, λ = h out at x = L .
This represents an approximation to a simple crushing switch such that rock fragments reduce in size given pressure larger than the threshold p crush , which is given below, and a lack of proximity to the inlet, where x crush defines the transition from the inlet zone to the crushing zone. Here, δ is our control parameter to adjust the rate of crushing, T c is a representative travel time through the crusher, and we take results from [35][36][37][38] in our derivation of the λ-dependent p crush given by We have assumed rough symmetry in the grains in defining p crush , little plasticity for the contact displacement s c and point out that the constant c is the modifier found in [37] to be commonly in the range As a result of the above equations, we have a system which features pressure dependence via these principal features: -A pressure-thickening material, -Effective pressure-roughening at the boundaries, -Pressure-dependent crushing.

Variational form
We now derive variational forms for use in the finite element approximation of the underlying system of partial differential equations. We defineü,v,p,λ as test functions associated with our unknowns. Concentrating first on our fluid problem, we multiply Eq. (15a) byp and (15b) by (ü,v), integrate over the channel Ω, and then apply integration by parts to obtain We use the surface integrals to prescribe our natural boundary conditions on τ and p, but must enforce zero flux across the physical boundaries via additional terms [39] in the variational formulation: where ξ is given by The new quantity q is the polynomial degree ofü andv, while h e is a piecewise constant function defined over Ω, which indicates the diameter of each element in the mesh. We choose approximating functions of the form: u, v piecewise bilinears plus cubic bubbles, p, λ piecewise bilinears, and select q = 2 as a compromise between the linear and cubic portions of our approximating functions and as a value in line with the ideal but computationally expensive biquadratic functions normally used to satisfy the stability condition. The variational forms for the fluid problem upon rewriting in terms of τ and p are therefore where Λ ∈ [0, 1] is a continuation parameter which we reduce during the calculation to achieve our normal stress boundary conditions at the inlet and outlet. Moving on to the particle size problem, we discretise Eq. (17) in time, by effectively approximating for small time step t. We therefore obtain Particle size λ is paired with its test functionλ from a piecewise bilinear finite element space and the system integrated over the channel Ω, to arrive at (21) We retain the initial and boundary conditions already stated.

Procedure
We compute our numerical solutions using Newton's method, and explicit time-stepping for the λ-problem. We employ sufficiently small time steps, deduced from typical values in Table 2, such that we will not have problems  with the Courant Friedrichs Lewy (CFL) condition. Specifically, we make an estimate such that material should not travel more than an average element diameterh e in distance over a time step t. We adopt the procedure illustrated in Fig. 15, noting that the initial continuation from the parallel-sided to nonparallel-sided channel problem includes continuation on five parameters simultaneously in a vector G: and represents not only continuation on geometry but also in the pressure drops down the channel and on normal stresses at the inlet and outlet. We prescribe an initial zero pressure drop along the channel and Λ = 1, and continue to finite pressure drop p step ≥ 0 as well as Λ = 0.

Results
We define the following functionals:  which is satisfactorily similar to the typical reference nondimensionalised values shown in Table 2. We establish a set of reference points down the channel, as illustrated in Fig. 16. Some plots from our reference solution are shown in Figs. 17, 18, 19, 20, 21 and 22, illustrating the steady pattern of that solution, which we believe to be numerically stable and representative of the physical expectations of our system. We see staggered crushing down the channel from the λ traces at the reference marks shown in Fig. 16 and the expected reduction in Q for increasing pressure and the same behaviour in the velocity field. The pressure conditions at the inlet and outlet are maintained at the level of the steady solution, and the flux across the slip boundary is appropriately minimised. Slippages at the basal and upper surfaces are seen to be successfully in antiphase with the applied stress. Upon viewing the pressure reference traces, there is a pressure hump mid-channel which is somewhat counter-intuitive, but could be explained as a necessary part of the solution to allow the required effective viscosities in the converging channel geometry. We now assess the impact of modifying the frequency and strength of the crushing stroke on the output of the crusher.

Varying stroke strength
We now investigate the effect of altered crushing stroke strength on the output quality and quantity. We employ continuation at t = 50 and calculate solutions for p h,high = 1.0, 1.1, 1.2, 1.3, 1.4. The bottom line results are plotted in Fig. 24 and seem to illustrate that increasing the crushing strength reduces not only the average grain size of the output but also the average speed of the material. The consistency measures both reduce with increasing p h,high , reflecting more homogeneous output material and less variance in the velocity field from its mean value over the time span.

Discussion
In this article, we have formulated a two-dimensional system and its variational equivalent, which approximates the flow of granular material in a two-dimensional fixed geometry cone crusher's cross section. By manipulating the model parameters as set out for the viscoplastic constitutive law and adding a simple crushing rule, we can approximate the flow of granular material and crushing under compressive crushing pressure in the channel. The crushing stroke in the fixed two-dimensional geometry is simulated by an oscillatory applied normal stress at the inlet and outlet, which can be modified in both frequency and size to assess the influence of those factors on the output size and consistency.
From our initial investigations, it seems that adjusting the frequency of the crushing stroke has very little effect on the quality or quantity of the output. There is a somewhat surprising result of improved consistency of output for small frequencies which should be investigated further, which is potentially linked to the fact that for such small frequencies, the stroke period approaches that of the transit time of the material in the crusher. Predictably, increasing the crushing strength does reduce the mean output size and improve the consistency of the output. This greater degree of crushing is offset by a reduced throughput of material.
In assessing the model, there are several outstanding issues. We have prescribed a pressure drop down the channel in order to force a sufficiently large throughput of material through the machine. For a realistic gravitydriven approximation to the real-life crusher situation, we need to investigate methods for facilitating such flows of granular media without the pressure difference. Increasing the inertial reference I 0 can suffice to a certain extent but saturates quickly. Reductions in the friction coefficient ζ w can also boost throughput, but very small ζ w and p h can cause computational problems due to infinite slippage and a model that requires non-negative pressure. The pressure 'hump' that exists mid-channel will also retard flow via the slip boundary conditions at the top and base as well as thickening the medium. The pressure satisfaction tests show errors on the order of 10 −6 at the inlet and outlet, which are to be considered satisfactory. A further simplification has been in the minimal choice of particle size dependence via the inertial number I , which should be re-evaluated in future research. Overall, this work provides a promising foundation for future work while being a ground-breaking first attempt to model the flow in a rock crusher using a rational, continuum approach.
Future work will be devoted to the three-dimensional generalisation of this problem, beginning with the fixed geometry problem and moving on to the moving boundary version.