A robust numerical method for granular hydrodynamics in three dimensions

Abstract Granular systems are of a discrete nature. Nevertheless, it can be advantageous to describe their dynamics by means of continuum mechanical methods. The numerical solution of the corresponding hydrodynamic equations is, however, difficult. Therefore, previous numerical simulations are typically geared towards highly specific systems and frequently restricted to two dimensions or mild driving conditions. Here, we present the first robust general simulation scheme for granular hydrodynamics in three dimensions which is not bound to the above limitations. The performance of the simulation scheme is demonstrated by means of three applications which have been proven as notoriously difficult for numerical hydrodynamic description. Although, by construction, our numerical method covers grain-inertia flows, the presented examples demonstrate that it produces reliable results even in the jammed or high density limit.

In three dimensions, the particle number density is related to volume fraction φ, and diameter of grains σ , via n = 6φ/(πσ 3 ). On Navier-Stokes level, their evolution obeys the hydrodynamic equations (Jenkins & Savage 1983) ∂n ∂t + ∇ · (nu) = 0, n ∂u ∂t + (u · ∇)u = − 1 m ∇ · P + ng, n ∂T ∂t + (u · ∇)T = −∇ · q − P : ∇u − ζ nT, where P is the pressure tensor, q is the heat flux, g is an external body force per unit mass, m is the mass of a single grain and ζ describes the rate of energy loss due to dissipative interaction of particles. To Navier-Stokes order, the components of the pressure tensor are given by P ij = pδ ij − τ ij , where p is the hydrostatic pressure and is the stress tensor. Here, η indicates the shear viscosity, γ the bulk viscosity, δ ij is the Kronecker symbol and A : B ≡ trace(A · B) = A ij B ji , where we apply the Einstein notation and sum over indices that appear twice. The granular heat flux is given by with the coefficient of thermal conductivity κ and where μ relates heat transfer to the gradient of density. Equations (2.1)-(2.3) show that the hydrodynamic simulation of granular systems involves two essential elements: first, a set of constitutive relations is needed which relates the hydrostatic pressure, p, the viscosities, γ , and η, the thermal coefficients, κ and μ, and the cooling coefficient, ζ , to the field variables and to material characteristics. Second, a numerical scheme to solve the Navier-Stokes equations (2.1) is needed, which is stable for all density and flow regimes a granulate can adopt. In the following, we will derive a numerical scheme which fulfils the above requirements.
Regarding the constitutive relations which are required to solve (2.1), we apply the constitutive relations by Garzó & Dufty (1999). These constitutive relations contain an Enskog factor to correct for the finite volume of granular particles. We use the finite volume correction derived by Torquato (1995). To our knowledge, this set of constitutive relations yields the most accurate description currently available. Other constitutive relations have been derived by Jenkins & Richman (1986) and Luding (2009), where the former holds for nearly static systems while the latter is only available for two-dimensional systems. In the dense limit close to random close packing there further exist incompressible approaches which assume constant volume fraction in the dense regime, thus, neglect small but important fluctuations of density due to irregular packings of granular materials, see e.g. Silbert, Landry & Grest (2003), Jop, Forterre & Pouliquen (2006) and Barker & Gray (2017). Moreover, the energy dissipation is not considered to be a consequence of binary dissipative particle collisions, thus, it is not compatible with kinetic theory. The latter point was improved by Barker & Gray (2017), who proposed a model that smoothly combines theories for dense and dilute flows. While these approaches may lead to excellent agreement with experiments and particle simulations, they do not provide information on the packing fraction or on velocity fluctuations. It has been shown that the incompressible μ(I) approaches suffer from ill posedness at high and low inertial numbers (Barker et al. 2015). Later, several well-posed viscous compressible theories were developed, see e.g. Heyman et al. (2017) and Schaeffer et al. (2019), as well as non-local theories that introduce further field variables to regularize the models, see e.g. Kamrin & Koval (2012) and Kamrin & Henann (2015).
It is advantageous to rewrite the hydrodynamic equations (2.1) in terms of the conservative variables density (n), momentum density (nu) and energy density (E = 1 2 mu 2 + 3 2 T) ∂n ∂t + ∇ · (nu) = 0, where I is the identity matrix and ⊗ is the tensor product. Albeit the sets (2.1) and (2.4) are mathematically equivalent, it can be shown (Hou & LeFloch 1994;Toro 2013) that the numerical solution of the non-conservative formulation (2.1) does not converge to the exact solution in the presence of shock waves. Therefore, it is inevitable that we work with the conservative formulation (2.4). The set of equations (2.4) resemble the corresponding set for a compressible, thermally conducting viscous fluid. It is, however, problematic to equate granular flows with compressible, thermally conducting viscous fluids due to the energy source term − 3 2 ηξ T. This term describes the inelastic interaction of the particles of the granulate and renders standard computational fluid dynamics (CFD) techniques inapplicable.
In order to rewrite the conservative formulation (2.4) in compact form, we introduce the five-dimensional vector of conservative variables, where the lower index of a vector indicates the corresponding vector element, e.g. u x stands for the x-component of u (note that only vectors with three spatial dimensions are printed in bold). We further introduce the convective fluxes y, z}, (2.5c) and the source term With these definitions, the conservative formulation (2.4) assumes then the compact form With this definition, ∇ · F is a 5 × 1 column vector like U and S(U).

Operator splitting
Equation (2.6) contains convective and diffusive contributions. In general, of course, systems reveal both convective and diffusive effects simultaneously, and this is the case for most granular systems of interest. Therefore, a universal numerical scheme must solve the convective and the diffusive equations simultaneously. The contributions are, however, of a very different nature: while diffusive processes affect the field variables according to their gradients in all directions, convection propagates the fields only in the flow direction (Versteeg & Malalasekera 2007). Therefore, the convective and diffusive contributions behave mathematically very differently and require very different numerical treatment, which causes difficulties when handled by the same numerical scheme (Holden et al. 2010): numerical schemes that are robust and efficient for the convective equations are frequently inappropriate for the diffusive contributions and vice versa. Therefore, we separate the hydrodynamic equations in such a way that convective and diffusive contributions can each be solved with individual, highly specialized numerical schemes. The convective contribution to the hydrodynamic equations (2.4) reads (3.1a,b) and the corresponding diffusive problem is given by Assume now an exact solution operator S c t to (3.1a,b) and an exact solution operator S ν step, t, such that To second order of accuracy, the solution of the combined convection-diffusion problem, (2.6), is then given by according to the splitting scheme by Strang (1968). Regarding the numerical solution of the hydrodynamic equations (2.4), the operators S c t and S ν t can be represented by specifically tailored numerical schemes which we will describe in the following sections. For a detailed introduction to operator splitting schemes, see Holden et al. (2010).

Solution operator for the convection equation
To derive an appropriate solution operator of the convective problem, (3.1a,b), we apply a finite volume method (FVM) where the simulation domain is decomposed into discrete volumes V. In contrast to frequently used finite difference discretization schemes, FVM inherently ensure the conservation laws. The surface of the volume elements is denoted by S. In practice we use the cells of the Cartesian mesh as volume elements and decompose the simulation domain into N (4.1) The surface of the volume elements is denoted by S. The key idea is now to transform a partial differential equation into a set of algebraic equations by integrating the partial differential equation over the finite volumes. The finite volume integration of (3.1a,b) reads Using the divergence theorem the volume integral can be replaced by an integral over the surface of the cell, where we additionally exchanged integration and differentiation and defined the cell average of the conservative variablesŪ ≡ 1/ V V dV U. The vector n denotes the unit vector normal to the surface, pointing outward. In the following we consider a cell indicated by C indicated by subscript e, w, etc. and variables referring to the body of the neighbouring cells by subscript E, W, etc. With these definitions, the surface integral in (4.3) expands to whereF c i,j with j ∈ s denoting the spatial average of the ith component of the convective flux over the cell face j, at time t. For the east cell face, we have, for instance, In the subsequent derivations we restrict ourselves toF c x,e , the remaining fluxes are treated analogously.
We approximate the integrals in (4.5) by Gaussian quadrature with two sampling points in each direction as illustrated in figure 2 where the arguments y α and z β are the coordinates of the integration points on the east surface and ω α , ω β are the corresponding weights. Equation (4.6) involves the values F c x,e of the convective flux at the quadrature points, which, in turn, depend on the conservative field variables U according to (2.5b). Indeed, we would need the values of U at the quadrature points, however, according to our finite volume discretization, (4.4), only the cell averages of U are available. Therefore, we interpolate the values of U at the quadrature points from cell averages of U of the cell where the quadrature points are located and the values of U in a certain neighbourhood of this cell. This interpolation is delicate as it determines the shock capturing abilities of the numerical scheme. Here we use a fifth-order weighted essentially non-oscillatory (WENO) scheme (Liu, Osher & Chan 1994) that provides an acceptable balance of accuracy and efficiency and delivers fifth-order accuracy in smooth regions of the field and third-order accuracy near discontinuities. For the one-dimensional case, WENO is explained in detail by Liu et al. (1994). For the extension to three-dimensional problems we additionally apply the method of dimension-by-dimension interpolation by Titarev & Toro (2004).
The quadrature points on the east face of a cell centred at x i coincide with the quadrature points on the west face of the neighbouring cell centred at x i+1 . However, the interpolation of U at the east face of a cell centred at is different from the interpolated value at the west face of the neighbouring cell at x i+1 (U R ≡ U w (x i+1 , y α , z β )), because a different stencil is considered for each interpolation. This deviation leads to a violation of conservation laws because U L / = U R implies that the flux leaving the cell centred at x i at its east face differs from the flux that enters the neighbouring cell centred at x i+1 at its west face. Therefore, the difference between U L and U R needs to be corrected. As shown e.g. in Toro (2013), the naïve approach of using the average of U L and U R leads to instabilities. Therefore, the flux between two neighbouring cells has to be derived from more elaborate upwind schemes that account for the direction of the flow. Moreover, to suppress spurious oscillations near discontinuities, the upwind scheme should assure monotonic flux. The multi-stage scheme by Titarev & Toro (2004) fulfils both criteria. The key idea of this scheme is to predict the flux from the initial data U L , U R and then to evolve U L and U R iteratively via the governing equation (3.1a,b). Here we use the first-order centred (FORCE) flux as a predictor where U We terminate the above iteration after k steps once the relative error of the predictor flux F FORCE is less than one per cent and assume F = F (k) FORCE . We can now compute the right-hand side of (4.4). The numerical solution of the convective equation (3.1a,b) corresponds to the solution of the ordinary differential equation (4.4). For this purpose we apply the third-order total variation diminishing Runge-Kutta method derived by Gottlieb & Shu (1998). The corresponding time iteration reads where the operator L corresponds to the right-hand side of (4.4).

Solution operator of the diffusion equation
We T and analogously for ∂ y U and ∂ z U. The diffusion matrices D ij are derived in Appendix A. These matrices are not restricted to the structured meshes used in this work, they can also be applied for e.g. triangular or tetrahedral meshes. We define further ∇U ≡ (∂ x U, ∂ y U, ∂ z U) T and to write the diffusive flux in the form can then be written as Similar as for the convection equation, we map this partial differential equation to a set of algebraic equations by means of the finite volume approach. The finite volume integration of (5.4) reads where we again used the divergence theorem, to replace the first integral on the right-hand side over the volume element V by an integral over its surface S. The diffusion matrices can be split into diagonal and non-diagonal (cross-diffusion) contributions. For the case of a Cartesian mesh, Darwish & Moukalled (2009) have shown that only the diagonal terms, D ii , contribute to the solution. In this case, the surface integral over the diffusive fluxes can be written as where S i is the surface element with the normal pointing in the i-direction, e.g. S x = y z.
If we expand the sum on the right-hand side over all faces and use a central difference approximation for the derivatives, the diffusive equation reads For convenience, we define and rewrite (5.7) in compact form The ordinary differential equation (ODE) in (5.9) is stiff in regions of high packing fractions, where the kinetic coefficients such as, e.g. the viscosity diverge. This is particularly relevant for driven systems where granulates frequently change in the course of time between very rarefied and random close packing states. The application of explicit methods is problematic for stiff ODEs as an impractically small time step may be required to keep the numerical error bound. Explicit methods would, thus, significantly restrict the range of applicability of the numerical scheme. Therefore, we use an implicit Runge-Kutta method (Butcher 1964) to solve (5.9) for all grid cells. In general, a s-stage Runge-Kutta scheme for the numerical solution of an initial value problem denotes the sampling points where the right-hand side of the ODE (5.10a,b) is evaluated. In general, (5.12) is a nonlinear system of equations which must be solved numerically to obtain the sampling points k i . To this end, we use the modified Newton iteration scheme by Houbak, Nøsett & Thomsen (1985). A scheme due to (5.12) can be specified by the coefficients a ij and b i which can be represented in a Butcher (1964) We implemented a two stage Gauss-Legendre quadrature formula which provides fourth-order accuracy. The coefficients of this scheme read The efficiency of the algorithm can be improved significantly by switching between this highly accurate scheme and a less accurate but more efficient scheme in regions where (5.9) does not reveal stiff behaviour. Such a switching between a highly accurate but numerically expensive algorithm and a numerically inexpensive implicit midpoint rule was proposed by Nørsett & Thomsen (1986). The implicit midpoint rule delivers second-order accuracy and its coefficients read in Butcher notation 1 2 1 2 (5.15)

Boundary conditions
Three types of boundary conditions are especially relevant for granular systems: periodic boundaries, solid walls and oscillating walls. While the implementation of periodic boundaries is straightforward, there are different ways to model solid boundaries. First, one can extend the simulation domain by ghost cells. The hydrodynamic variables of these ghost cells are determined from the corresponding values of their neighbouring cells within the physical simulation domain whilst taking into account the mirror symmetry with respect to the bounding wall (Carrillo et al. 2008). For instance, to model a slip boundary condition with no mass flux through the wall, the normal component of the velocity field is inverted while the tangential component of the velocity field, as well as the energy field and the density field, are copied. Similarly, a no-slip boundary can be achieved by additionally inverting the tangential component of the velocity field. Another way to model boundary conditions is to determine the hydrodynamic variables on the boundary explicitly (Blazek 2015). This approach reduces the numerical errors due to the interpolation of the fluxes on the boundary surface. It is, thus, more consistent with the applied finite volume method. To represent slip boundaries, only the normal component of the velocity field is set to zero at the boundary, while for a no-slip boundary all components of the velocity field are set to zero at the boundary. For the no-slip solid wall the convective and diffusive fluxes in (2.5b) and (2.5c), reduce to where p w denotes the pressure at the wall which can be extrapolated from the neighbouring cells within the simulation domain (Blazek 2015) and q i is the heat flux through the boundary. In case of elastic walls, there is no energy loss and the heat flux q i vanishes at the boundary. If the interaction between the boundary and the particles of the granulate is dissipative, there is a finite heat flux through the boundary. This heat flux is determined by the mean collision frequency, ν 0 , at which the particles collide with the boundary (Brilliantov & Pöschel 2010) and by the coefficient of restitution ε w which describes the fraction of energy that is dissipated when a grain collides with the wall where A s denotes the surface area of the boundary. Equation (6.2) is obtained by multiplying the mean energy loss of a particle impacting the wall by the mean collision frequency per surface area. In many cases, granular systems are driven by mechanically vibrating walls with a certain frequency and amplitude. The corresponding oscillating walls can be modelled by rewriting the hydrodynamic equations (2.4) in the non-inertial frame of the vibrating walls (Carrillo et al. 2008). The corresponding inertial forces act on all particles in the system, irrespective of their position. Therefore, in case of driven systems, inertial forces have to be added to the momentum equation in (2.4).
The interaction between granular matter and the system boundaries can be much more complicated. The grains can, e.g. roll along the boundary surface. Additionally, the shape and the surface roughness of the grains have great impact on the boundary-bulk interaction. However, these grain-level details are not subject of a hydrodynamic model.

Numerical validation
Because of the complexity of the hydrodynamic equations, (2.4), a theoretical validation of the presented numerical scheme is not feasible. To our knowledge, there are no experimental results, in which temperature, velocity and density are measured simultaneously, which precludes us from an experimental validation of our approach as well. Therefore, we validate our hydrodynamic simulations against well-established particle-based simulations of the corresponding granular system. For the validation we consider the following scenarios: (i) a granular gas which sediments under the action of gravity; (ii) a bed of grains in a vertically shaken container in presence of gravity; and (iii) the flow of granular material down an inclined pile.
These set-ups are demanding, regarding a hydrodynamic simulation because they involve extreme variations of the packing fraction, from dilute gas to random close packing, and steep gradients in the hydrodynamic fields. For the particle-based reference simulations we apply event driven molecular dynamics (EDMD) and the discrete element method (DEM) (see, e.g. Pöschel & Schwager (2005), and many references therein). Both methods describe the granular systems reliably on different levels of detail. An EDMD simulation is considered as the reference for this test. In EDMD simulations, the particles are modelled ideally rigid. A collision between two ideally rigid, smooth particles corresponds to an instantaneous transfer of momentum, that is, to an instantaneous change of the particles' velocities. Collisions are characterized by the coefficient of restitution, ε relating the pre-and post-collisional values of the normal components of the particles' relative velocities at the point of contact. In our simulation we assume ε = 0.9. Further details on event-driven molecular dynamics can be found e.g. in Bannerman, Sargant & Lue (2011) and references therein.

Sedimentation test
To compare the results of the hydrodynamic simulations with the results of the particle-based EDMD simulations, continuous hydrodynamic fields, such as density, flow velocity and energy density, have been computed from the discrete particle data. For the sedimentation test, the system is homogeneous in the horizontal direction while gravity acts in the vertical direction, therefore, we can represent the fields as functions of the vertical component (height). Figure 3 compares the conservative variables (density, momentum density and energy density) as functions of the height (direction of gravity) as obtained from the hydrodynamic simulations with those obtained from the reference EDMD simulations. We see that the field of density is well described by the hydrodynamic approach for all times and heights. For the fields of velocity and energy, the hydrodynamic description and the particle-based reference simulation are in good agreement for times up to approximately 100 ms. For later times, we still see qualitative agreement of the curves, however, the values deviate moderately. Similar behaviour can be observed for the field of temperature shown in figure 4. The numerical results are not affected by substantial refinements of the spatial and temporal resolutions of the numerical scheme which indicates the reliability of our numerical solution. We therefore attribute the deviations to deficiencies of the underlying continuum mechanical model. The constitutive relations provided by Garzó & Dufty (1999) and the applied equation of state are approximations to the physics of granulates. As the applied constitutive relations are derived from kinetic theory, they are only valid for dynamic systems. The deviations start to get significant when part of the grains rest on the floor (t = 150 ms). Additionally, the hydrodynamic equations we solve are approximated to first-order derivatives which may cause error when sharp gradients in the hydrodynamic fields are present. To improve the accuracy, one should apply hydrodynamic equations beyond the Navier-Stokes level and incorporate terms of Burnett order and super-Burnett order (Chapman & Cowling 1939;Sela, Goldhirsch & Noskowicz 1996;Sela & Goldhirsch 1998;Goldhirsch 1999Goldhirsch , 2001  Despite the fully three-dimensional simulations which have been applied, the present sedimentation test is effectively a one-dimensional problem. However, to our knowledge, it is the first test case in which all relevant variables (density, momentum and energy) as obtained from a hydrodynamic simulation are compared with the corresponding results from the particle system.

Bouncing bed test
We consider granular particles under the action of gravity in a vertically shaken container. This system exhibits various modes of behaviour. In an elaborate experimental study, this rich behaviour was compiled to a phase diagram by Eshuis et al. (2007). For a quantitative evaluation of our hydrodynamic simulations, we compare the hydrodynamical results with those obtained from corresponding particle-based discrete element simulations, that is, we solve Newton's equation of motion for the translational and rotational motion of all N grains.
When particles collide, they mutually deform each other, quantified by and a tangential component f t perpendicular to e n f = f n + f t (7.2) For viscoelastic particles, the normal component reads (Brilliantov et al. 1996) where Y and ν denote Young's modulus and Poisson's ratio, respectively, and R eff = R i R j /(R i + R j ). The dissipative parameter A dis is a function of the material viscosity, see Brilliantov et al. (1996). As the coefficient of restitution is a (known) function of the interaction force (Schwager & Pöschel 1998a,b;Ramírez et al. 1999) and the interaction velocity, the coefficient A dis can be determined provided the coefficient of restitution for a particular collision rate is known (Müller & Pöschel 2011), for instance from an experiment.
The tangential component of the interaction force causes a torque on the colliding particles, which, in turn, alters their angular velocity. Here, we apply the model by Cundall & Strack (1979) which mimics static and dynamic friction by means of an imaginary spring (for an extended discussion see Pöschel & Schwager 2005;Schwager, Becker & Pöschel 2008). The spring of initial length zero is created at time t c when the particles touch. Its elongation is then where v t rel ≡ e t · [ṙ ij + e n × (R i ω i + R j ω j )] is the tangential component of the relative velocity of the particles at the point of contact and ω i and ω j are the angular velocities of particles i and j. The tangential component of the interaction force is then given by the restoring force of the spring, limited by Coulomb's law f t = − sign(v t rel ) · min(|KZ| , μ| f n |). (7.8) Here, μ denotes Coulomb's coefficient of friction and K is the elastic constant of the spring. In contrast to the normal force, the parameters of the tangential force, K and μ, cannot be determined from material bulk properties and geometry but they are mainly determined by surface properties. In practice, they can be only determined by comparison of the results of benchmark simulations and experimental results. For the bouncing bed test, we consider a rectangular box with of size L x × L y × L z = 100 × 100 × 10 mm 3 , filled by 8560 particles of diameter 1 mm. The material parameters used in the simulation read Y = 5 × 10 9 Pa, material density 2.6 g cm −3 , ν = 0.3, A dis = 1.61 × 10 −7 s, K = 0.347ρ el R eff and μ = 0.7. The particles are vertically vibrated in the presence of gravity with the frequency 12 Hz and the amplitude 4 mm. For the interaction between particles and the wall we assume the same parameters. The integration time step is 10 −5 s. To compare the results of the DEM simulation with the results of the corresponding hydrodynamic simulations, the domain is divided into 20 × 100 × 4 Cartesian grid and a time step of 1/12 × 10 −5 s is used. We consider the time-dependent density field as a function of the vertical coordinate, again averaged over the homogeneous horizontal dimensions. We assume no-flux boundary conditions at the upper side of the simulation domain near the free surface. Figure 5 shows the packing fraction as a function of the distance from the bottom wall as obtained from the hydrodynamic simulations in comparison with the results from the corresponding DEM simulation. The results of both methods are in good quantitative agreement for all phases of the vibration cycle. Similar results were obtained before by e.g. Bougie, Policht & Pearce (2012). There, however, mild driving conditions are considered which lead to maximum volume fractions way below 0.3. In our study we apply driving parameters causing packing fractions larger than 0.5. Due to the extreme conditions of the bouncing bed set-up (dynamic transitions between densely packed and dilute regions) the temperature cannot be reliably obtained from particle simulations. For high density where the granulate occurs in an almost crystalline or random close packed state, the relative velocities of the grains obtained from particle simulations are dominated by unavoidable numerical noise, which impedes a meaningful evaluation of the granular temperature. For the dilute extreme, the number of particles in these regions is almost zero; see, e.g. figure 5(b) for small height. Thus, in extremely dilute regions, we cannot compute the value of granular temperature from particle simulation data. Nevertheless, the temperature as obtained from granular hydrodynamics can be computed also in the dilute and the densely packed regions, and this is shown in the bottom raw of figure 5.

Flow on a pile with rough sidewalls
The development of general constitutive relations for the continuum description of granular flows is still work in progress. The constitutive relations for the gas regime can be derived from kinetic theory (see, e.g. Goldhirsch 2003) and for the solid regime from soil mechanics (see, e.g. Nedderman 1992). A constitutive law for dense granular flow is proposed by Jop et al. (2006), where the granular flow is described as incompressible with a simple visco-plastic relation for viscosity. The law of Jop et al. (2006) is rate independent to the extent that only the coefficient of friction depends on the shear rate. While the constitutive law by Jop et al. (2006) leads to excellent agreement with experiments, it provides no information about packing fraction and velocity fluctuations inside the bulk.
In this section, we consider glass beads of diameter σ = 1 mm, in a channel of length 226 cm and width w = 48σ , with rough sidewalls. The thickness of the flowing material, h, is initialized to 22.5 mm. The channel is tilted by angle θ = 27 • . Figure 6 shows a sketch of the set-up. To compare our results with the simulations by Silbert et al. (2003) and Jop et al. (2006), we define the dimensionless flow rate Q = Q/(σ 3/2 g 1/2 ), where Q is the 917 A33-17  flow rate per unit of width w, g is the gravitational acceleration and V = V/(gσ ) 1/2 is the dimensionless velocity in the x direction.
The model by Jop et al. (2006) describes dense granular flow as an incompressible fluid (see discussion in § 2). Here, we solve the full Navier-Stokes equations (2.4) using the constitutive relations from Garzó & Dufty (1999). We discretize the domain into a 80 × 80 × 40 structured grid and use the time step 10 −6 s for the time integration. The rough sidewalls and the channel floor are considered as no-slip boundaries, as described in § 6. We assume no-flux boundary conditions at the upper side of the simulation domain near the free surface.
The first consequence of solving the full equations (2.4) is that the volume fraction is no longer constant. Figure 7 shows the profile of the packing fraction along the y direction.  (2001). We find good agreement between the experimental and numerical density profiles. The different asymptotic densities (Φ for small y) may be attributed to the fact that the experiment uses slightly polydisperse grains, whose density at random close packing is larger than for monodisperse spheres.) As expected, in the lower part of the flow, the material is more dense, while there is a dilute gaseous region on the surface of the pile. The existence of this gaseous region is in good agreement with the experiment performed by Forterre & Pouliquen (see, e.g. figure  1 in Forterre & Pouliquen 2008). For a more quantitative comparison, figure 7 shows experimental results by Ancey (2001) for channel flow of glass beads of 1 mm diameter down a rough incline (width 48 mm, slope 27 • , layer thickness 26 mm), see figure 17 in Ancey (2001). We find good agreement between the experimental and numerical density profiles. The different asymptotic density (Φ for small y) may be attributed to the fact that the experiment uses slightly polydisperse grains ('poorly sorted', no interval given, see table 1 in Ancey 2001). Note that the width of the numerically obtained profile in figure 7 extends to y ≈ 26 mm, which is larger than the initial layer width, h, given above. This inflation can be explained by Reynolds dilatancy of the flowing sand. The initial height, h, has been chosen such that the width of the stationary flow reaches approximately the same extension as the experimental data. The three-dimensional velocity profile resulting from the hydrodynamic simulation is shown in figure 8. Due to the rough sidewalls and the rough base of the pile, the velocity is lower close to these surfaces. In our simulation, the velocity profile decays slower as compared to the simulations by Jop et al. (2006). The reason, therefore, is that the viscosity is independent of pressure and only depends on temperature and filling fraction in our constitutive relations. On the contrary, the viscosity depends on pressure in the constitutive relations by Jop et al. (2006). As the pressure is higher at the bottom of the channel, the viscosity diverges to infinity faster. As a result, lower velocities are observed close to the bottom of the channel in Jop et al. (2006).

Summary
We presented a numerical scheme for the hydrodynamic simulation of granular systems in three dimensions. To our knowledge, this is the first general three-dimensional simulation algorithm for granular hydrodynamics which is neither restricted to certain dynamic states of the granulate nor to special boundary conditions. Most importantly, unlike previous methods used, e.g. in soil mechanics or rapid gas flows, our method is not restricted to a certain range of density. While, by construction, our numerical method covers grain-inertia flows, it was demonstrated that the algorithm is numerically robust for all values of packing fraction, from dilute gas to random close packing and likewise also with respect to harsh driving. However, one has to distinguish between the above-mentioned capabilities of the numerical scheme used to solve the Navier-Stokes equations, on the one hand, and the validity of the applied constitutive relations, on the other. The fundamental prerequisites to successfully apply the presented algorithm are, of course, constitutive relations which are valid over the whole range of densities and flow regimes. While massive advances have been made in this direction over the past twenty years, this ultimate research objective is still unreached. Comparison of the presented numerical scheme against well-established particle-based simulation techniques reveals two important points: (i) There are situations where the hydrodynamic results are in very good qualitative and quantitative agreement with the particle-based reference simulations. In these cases, both simulation methods support one another and we can assume that both methods are correct, including the underlying hydrodynamic model and the numerical scheme presented here. (ii) There are situations where we find qualitative agreement between the hydrodynamic and particle simulations, however, non-negligible quantitative deviations between the numerical data gained by hydrodynamics and particle mechanics. This, in turn, may be due to three different arguments: (a) The underlying hydrodynamic model is inadequate, that is, either the constitutive equations or the equation of state are inappropriate for the investigated system. This is not surprising as all available constitutive relations and equations of state are approximative and based on phenomenological arguments rather than on rigorous derivations. (b) Particularly for systems at large packing density, the Navier-Stokes level of the hydrodynamic description may be insufficient, that is, linear relations for the transport coefficients may be insufficient for the description of granular flow at high density. In this case, higher-order descriptions on Burnett or super-Burnett order must be considered, see the discussion in § 7.1. (c) The hydrodynamic description itself may be inadequate due to the arguments related to the lack of scale separation discussed at the end of § 1. In these cases, there is no hope for a hydrodynamic description. Fortunately, sufficiently dynamic systems of practical interest seem not to belong to this class. Our demonstration examples have been chosen to be hard for a hydrodynamic description, but still we could obtain stable results in approximate quantitative agreement with particle simulations.
It would be desirable to compare our results with competing approaches to solve the equations (2.4). However, all available studies apply explicit schemes and fail in the dense regime. A comparison which is restricted to the dilute limit is, however, inappropriate as the aim of the presented implicit method is neither to increase performance nor to optimize accuracy but only to improve the stability of the solver and, finally, to enable simulations in the dense regime.
The computational granular hydrodynamics presented here can be extended and generalized to other constitutive relations (e.g. Jenkins & Richman 1986), and equations of state (e.g. Luding 2009), available in the literature. The reliable and numerically stable hydrodynamic simulation method finally allows us to study large-scale granular systems which are not accessible through particle simulations due to the large number of particles involved, such as sand dunes (e.g. Schwämmle & Herrmann 2003), planetary rings (e.g. Grenberg & Brahic 1984;Esposito 2006;Fridman & Gorkavyi 2010), avalanches (e.g. Pudasaini & Hutter 2007) or large-scale technical apparatus (e.g. Savage & Hutter 1989).