Underwater granular flows down inclined planes

In this study, two-dimensional sub-grain scale numerical simulations are performed to understand the local rheology of dense granular flows in fluid. The Discrete Element (DEM) technique is coupled with the Lattice Boltzmann Method (LBM), for fluid-grain interactions, to understand the evolution of immersed granular flows. The fluid phase is simulated using Multiple-Relaxation-Time LBM (LBM-MRT) for numerical stability. The Eulerian nature of the LBM formulation, together with the common explicit time step scheme of both LBM and DEM makes this coupling strategy an efficient numerical procedure for systems dominated by both grain-grain and grain-fluid interactions. In order to simulate interconnected pore space in 2D, a reduction in the radius of the grains (hydrodynamic radius) is assumed during LBM computations. By varying the hydrodynamic radius of the grains, granular materials of different permeabilities can be simulated.A parametric analysis is performed to assess the influence of the granular characteristics (initial packing, permeability, slope of the inclined plane) on the evolution of flow and run-out distances. The effect of hydrodynamic forces and hydroplaning on the run-out evolution is analysed by comparing the mechanism of energy dissipation and flow evolution in dry and immersed granular flows.Voronoi tesselation was used to study the evolution of local density and water entrainment at the flow front. © 2015 Taylor & Francis Group.


INTRODUCTION
The flow of dense granular material is a common phenomenon in engineering predictions, such as avalanches, landslides, and debris-flow modelling. Despite the huge amount of research that has gone into describing the behaviour of granular flows, a constitutive equation that describes the overall behaviour of a flowing granular material is still lacking. The initiation and propagation of submarine granular flows depend mainly on the slope, density, and quantity of the material destabilised. Although certain macroscopic models are able to capture the simple mechanical behaviours, the complex physical mechanisms that occur at the grain scale, such as thydrodynamic instabilities, the formation of clusters, collapse, and transport, have largely been ignored (Topin et al. 2011). The momentum transfer between the discrete and the continuous phases significantly affects the dynamics of the flow (Peker and Helvac 2007). Grainscale description of the granular material enriches the macro-scale variables, which poorly account for the local rheology of the materials. In order to describe the mechanism of saturated and/or immersed granular flows, it is important to consider both the dynamics of the solid phase and the role of the ambient fluid (Denlinger and Iverson 2001). In particular, when the solid phase reaches a high volume fraction, it is important to consider the strong heterogeneity arising from the contact forces between the grains, the drag interactions which counteract the movement of the grains, and the hydrodynamic forces that reduce the weight of the solids inducing a transition from dense compacted to a dense suspended flow (Meruane et al. 2010). The case of the collapse in presence of an interstitial fluid has been less studied. In this paper, we study the submarine granular flows in the inclined configuration. We study the effect of permeability, density and slope angle on the run-out evolution.

LBM FORMULATION
The Lattice Boltzmann Method is a 'micro-particle' based numerical time-stepping procedure for the solution of incompressible fluid flows. Consider a 2D incompressible fluid flow with density ρ and kinematic viscosity v, in a rectangular domain D. The fluid domain is divided into a rectangular grid or lattice, with the same spacing 'h' in both the xand the ydirections, as shown in figure 1. The present study focuses on two-dimensional problems, hence the D2Q9 momentum discretisation is adopted (see He et al. (1997) for naming convention). Figure 1: The Lattice Boltzmann discretisation and D2Q9 scheme: (a) a standard LB lattice and histogram views of the discrete single particle distribution function/direction-specific den- The lattice Boltzmann Bhatnagar-Gross-Krook (LGBK) method is capable of simulating various hydrodynamics (Succi 2001) and offers intrinsic parallelism. Although LBM is successful in modelling complex fluid systems, such as multiphase flows and suspensions in fluid, the LBM may lead to numerical instability when the dimensionless relaxation time τ is close to 0.5. The Multi-Relaxation Time Lattice Boltzmann Method (LBM-MRT) overcomes the deficiencies of linearlised single relaxation LBM-BGK, such as fixed Prandtl number (Pr=ν/κ), where the thermal conductivity 'κ' is unity (Liu et al. 2003). The LB-MRT model offers better numerical stability and has more degrees of freedom. In the formulation of the linear Boltzmann equation with multiple relaxation time approximation, the lattice Boltzmann equation is written as: where S is collision matrix. The nine eigen values of S are all between 0 and 2 so as to maintain linear stability and the separation of scales, which means that the relaxation times of non-conserved quantities are much faster than the hydrodynamic time scales. The LGBK model is the special case in which the nine relaxation times are all equal and the collision matrix S = 1 τ I, where I is the identity matrix. The evolutionary progress involves two steps, advection and flux. The advection can be mapped to the momentum space by multiplying through by a transformation matrix M and the flux is still finished in the velocity space. The evolutionary equation of the multi-relaxation time lattice Boltzmann equation is written as: where M is the transformation matrix mapping a vector f in the discrete velocity space The transformation matrix M can be constructed via Gram-Schmidt orthgonalisation procedure. Through the Chapman-Enskog expansion (Du et al. 2006), the incompressible Navier-Stokes equation can be recovered and the viscosity is given as:

Turbulence in Lattice Boltzmann Method
Modelling fluids with low viscosity like water remains a challenge, necessitating very small values of h, and/or τ very close to 0.5 (He et al. 1997). Turbulent flows are characterised by the occurrence of eddies with multiple scales in space, time and energy. In this study, the Large Eddy Simulation (LES) is adopted to solve for turbulent flow problems. The separation of scales is achieved by filtering of the Navier-Stokes equations, from which the resolved scales are directly obtained and unresolved scales are modelled by a one-parameter Smagorinski sub-grid methodology, which assumes that the Reynold's stress tensor is dependent only on the local strain rate (Smagorinsky 1963). The turbulent viscosity ν is related to the strain rate S ij and a filtered length scale 'h' as follows: where S c is the Smagorinski constant found to be close to 0.03 (Yu et al. 2005).
The effect of the unresolved scale motion is taken into account by introducing an effective collision relaxation time scale τ t , so that the total relaxation time τ * is written as: where τ and τ t are respectively the standard relaxation times corresponding to the true fluid viscosity v and the turbulence viscosity v t , defined by a sub-grid turbulence model. The new viscosity v * corresponding to τ * is defined as: The Smagorinski model is easy to implement and the Lattice Boltzmann formulation remains unchanged, except for the use of a new turbulence-related viscosity τ * . The component s 1 of the collision matrix becomes s 1 = 1 τ +τt .

COUPLED LB -DEM FOR FLUID-PARTICLE INTERACTIONS
The Lattice Boltzmann approach has the advantage of accommodating large particle sizes and the interaction between the fluid and the moving particles can be modelled through relatively simple fluid -particle interface treatments. Further, employing the Discrete Element Method (DE) to account for the particle/particle interaction naturally leads to a combined LB -DEM solution procedure. The Eulerian nature of the Lattice Boltzmann formulation, together with the common explicit time step scheme of both the Lattice Boltzmann and the Discrete Element makes this coupling strategy an efficient numerical procedure for the simulation of particle-fluid systems (Cook et al. 2004). In order to capture the actual physical behaviour of the fluid-particle system, the boundary condition between the fluid and the particle is modelled as a non-slip boundary condition, i.e. the fluid near the particle should have similar velocity as the particle boundary. The solid particles inside the fluid are represented by lattice nodes. The discrete nature of lattice will result in stepwise representation of the surfaces. Very small lattice spacing is adopted to obtain smoother boundaries.

UNDERWATER GRANULAR FLOWS
In this study, a 2D poly-disperse system (d max /d min = 1.8) of circular discs in fluid was used to understand the behaviour of granular flows on inclined planes (see Figure 2). The soil column was modelled using 1000 discs of density 2650 kg m −3 and a contact friction angle of 26 • . The collapse of the column was simulated inside a fluid with a density of 1000 kg m −3 and a kinematic viscosity of 1 × 10 −6 m 2 s −1 . The choice of a 2D geometry has the advantage of cheaper computational effort than a 3D case, making it feasible to simulate very large systems. A granular column of aspect ratio 'a' of 0.8 was used. A hydrodynamic radius r = 0.9R was adopted during the LBM computations. Dry analyses were also performed to study the effect of hydrodynamic forces on the run-out distance.

Soil
Water g gate Figure 2: Underwater granular collapse set-up

Effect of initial density
The morphology of the granular deposits in fluid is shown to be mainly controlled by the initial volume fraction of the granular mass and not by the aspect ratio of the column (Rondon et al. 2011, Pailha et al. 2008. In order to understand the influence of the initial packing density on the run-out behaviour, a dense sand column (initial packing density, Φ = 83%) and a loose sand column (Φ = 79%) were used. The granular columns collapse and flow down slopes of varying inclinations (2.5 • , 5 • and 7.5 • ). The evolution of run-out distances for a dense sand column with time in dry and submerged conditions for varying slope inclinations are presented in figure 3. The run-out distance is longer in submerged condition than the dry condition for a flow on a horizontal surface. However, with increase in the slope angle the run-out in the fluid decreases.
Dense granular columns in fluid take a longer time to collapse and flow, due to the development of large negative pore-pressures, as the dense granular material dilates during the initial phase of the flow. The morphology of dense granular flows down slopes of varying inclinations at the critical time (t = τ c = H/g, when the flow is fully mobilised) are shown in figure 5.
It can be seen that the viscous drag on the dense column tend to predominate over the influence of hydroplaning on the run-out behaviour. This influence can be observed in the smaller peak kinetic energy for granular column in fluid compared to it's dry counterpart (see Figure 4). With increase in slope angle, the volume of material that dilates increases. This results in large negative pore pressures and more viscous drag on the granular material. Hence, the difference in the run-out between the dry and the submerged condition, for a dense granular assembly, increases with increase in the slope angle.
In contrast to the dense granular columns, the loose granular columns (relative density I D = 30%) show longer run-out distance in immersed conditions (see Figure 6). The run-out distance in fluid increases with increase in the slope angle compared to the dry cases. Loose granular material tends to entrain more water at the base of the flow front, creating a lubricating surface, which causes longer run-out distance (see Figure 7). The hydroplaning effect causes an increase in the velocity the loose condition in comparison with the dense condition (see Figure 8). The evolution of packing density (see Figure 9) shows that dense and the loose conditions reach similar packing density. This indicates that the dense granular column dilates more and is susceptible to higher viscous drag forces. Where as in the loose condition, a positive pore-pressure is observed at the base of the

Effect of permeability
In DEM, the grain -grain interaction is described based on the overlap between the grains at the contact surface. In a 3D granular assembly, the pore spaces between grains are interconnected. Whereas in a 2-D assembly, the grains are in contact with each other that result in a non-interconnected pore-fluid space. This causes a no flow condition in a 2-D case. In order to overcome this difficulty, a reduction in radius is assumed only during the LBM computation phase (fluid and fluid -solid interaction). The reduction in radius allows interconnected pore space through which the surrounding fluid can flow. This technique has no effect on the grain -grain interactions computed using DEM. See Kumar et al. (2012) for more details about the relationship between reduction in radius and permeability of the granular assembly. For a slope angle of 5 • , the hydrodynamic radius of the loosely packed grains was varied from r = 0.7R (high permeability), 0.75R, 0.8R, 0.85R to 0.9R (low permeability). The run-out distance is found to increase with decrease in the permeability of the granular assembly (see Figure 10). The run-out distance for high permeable conditions (r = 0.7R -0.8R) were lower than their dry counterparts. Although, decrease in permeability resulted in an increase in the runout distance, no significant change in the run-out behaviour was observed for a hydrodynamic radius of up to 0.8R. With further decrease in permeability (r = 0.85R and 0.9R), the run-out distance in the fluid was greater than the run-out observed in the dry condition. At very low permeability (r = 0.9R), granular material started to entrain more water at the base, which causes a reduction in the effective stress accompanied by a lubrication effect on the flowing granular media. This can be seen as a significant increase in the peak kinetic energy and the duration of the peak energy, in comparison with dry and high permeable conditions (see Figure 12).
The permeability of the granular column did not have an influence on the evolution of height during the flow. However, dry granular column tends to collapse more than the immersed granular column (see Figure 11).
Positive pore-pressure generation at the base of the flow was observed for low permeable conditions. Inspection of the local packing density showed entrainment of water at the base of the flow, which can also be observed by the steep decrease in the packing density (see Figure 13) for the very low permeability con- dition (r = 0.9R). At the end of the flow (t ≥ 3 × τ c ), the excess pore-pressure dissipates and the granular material, irrespective of their permeability, reaches almost the same packing density.

SUMMARY
Two-dimensional LB-DEM simulations were performed to understand the behaviour of submarine granular flows. Unlike dry granular collapse, the runout behaviour in fluid is dictated by the initial volume fraction. Granular columns with loose packing tend to flow longer in comparison to dense columns, due to entrainment of water at the base resulting in lubrication. The loose column when it starts flowing expands and ejects liquid, leading to a partial fluidization of the material. However, with increase in the slope angle, the run-out in fluid is influenced by the viscous drag on the granular materials. The run-out distance in fluid increases with decrease in permeability. More research work is required to characterise the flow behaviour of granular materials, especially in submerged conditions.