# Numerical Modeling of a Granular Collapse Immersed in a Viscous Fluid

Edouard Izard, Laurent Lacaze, Thomas Bonometti, Annaïg Pedrono

► To cite this version:

Edouard Izard, Laurent Lacaze, Thomas Bonometti, Annaïg Pedrono. Numerical Modeling of a Granular Collapse Immersed in a Viscous Fluid. Springer. Advances in Hydroinformatics SimHydro 2017: Choosing The Right Model in Applied Hydraulics, pp.1099 - 1116, 2018, 978-981-10-7218-5. 10.1007/978-981-10-7218-5_76. hal-04354015

HAL Id: hal-04354015

https://hal.science/hal-04354015v1

Submitted on 19 Dec 2023

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
**Official URL:**

https://doi.org/10.1007/978-981-10-7218-5_76

**To cite this version:**

Izard, Edouard and Lacaze, Laurent and Bonometti, Thomas and Pedrono, Annaig Numerical Modeling of a Granular Collapse Immersed in a Viscous Fluid. (2018) In: SimHydro 2017, 14 June 2017 - 16 June 2017 (Sophia Antipolis, France).
Numerical Modeling of a Granular
Collapse Immersed in a Viscous Fluid

Edouard Izard, Laurent Lacaze, Thomas Bonometti
and Annaïg Pedrono

# 1 Introduction

Dense particulate flows in an incompressible Newtonian fluid are encountered in a
large number of applications such as civil and chemical engineering, transportation,
steel and food industries, geophysics and oceanography. In these applications,
granular flows can be gravity-driven, shear-driven by a fluid flow or both. Here, we
focus on a gravity-driven case, namely the collapse of granular column immersed in
a Newtonian fluid where grains motion is induced by the difference of density
between the granular media and the surrounding fluid. Note that depending on the
fluid properties (density, viscosity), the characteristics of the flow can be drastically
different due to the relative influence of particles contact forces and fluid-particle
interaction forces. Moreover, the processes involved in these granular avalanches
from the particle scale to that of the full system are not fully understood. In
particular, it is complicated to experimentally measure fluid and particles local
dynamics inside a moving granular bed.

Empirical observations of granular columns collapse in the air were conducted
either in an axisymmetric configuration by Lube et al. [1] and Lajeunesse et al. [2]
or in two dimensions in a rectangular channel by Lube et al. [3] and Balmforth and

E. Izard (✉)

ArcelorMittal R&D Maizières, Voie Romaine, 57283 Maizières-Lès-Metz, France

e-mail: izard.edouard@gmail.com

L. Lacaze · T. Bonometti · A. Pedrono

IMFT, Université de Toulouse, CNRS, INPT, UPS, Toulouse, France

e-mail: Laurent.Lacaze@imft.fr

T. Bonometti

e-mail: thomas.bonometti@imft.fr

A. Pedrono

e-mail: annaig.pedrono@imft.fr

https://doi.org/10.1007/978-981-10-7218-5_76
Kerswell [4]. In both cases, the final granular deposit is shown to be a function of the initial geometry of the column, namely, the initial aspect ratio *a* defined as the ratio between the initial column height *H*<sub>o</sub> and its radius *R*<sub>o</sub>. Thompson and Huppert [5] have observed a strong modification in the collapse dynamics of a granular column when the ambient fluid play a role, e.g. using water as surrounding fluid. More recently, an experimental investigation of the granular collapse in a viscous fluid highlighted the strong influence of the initial packing fraction of the granular column on the granular dynamics and deposits [6]. In particular, the fluid pressure, relative to the hydrostatic pressure, measured at the basis of the column was shown to be positive (resp. negative) when the initial packing fraction of the column is below (resp. above) the critical value of Φ<sub>ic</sub> ≈ 0.59, at the initial stage of the collapse. In addition, this difference of initial dynamics alters the granular deposit at small aspect ratio between a dense and loose initial state of the granular media. For instance, for *a* = 0.5, the shape of the deposit can be trapezoidal or triangular in a dense or loose configuration, respectively. This mechanism, partly due to the contracting nature or dilatancy of the granular materials, is called pore pressure feedback by Iverson [7].

Numerical simulations of dry granular collapses, e.g. when the ambient fluid has no effect on the granular dynamics, using the discrete element methods (DEM) qualitatively reproduce the experimental observations (see [8–11]). Simulations of granular collapses, where the fluid has an influence on the grains motion, are rare. Topin et al. [12] performed two-dimensional simulations of granular collapses in a viscous fluid using a finite element method for the resolution of Navier–Stokes equations coupled with the non-smooth contact dynamics method for the grains motion. In this method, an artificial ring of fluid is prescribed around the particles so that the fluid can flow around the grains in a dense two-dimensional granular configuration. Fixing the initial packing fraction of the column, the authors observed deposition laws in the viscous, inertial and dry regimes. These regimes were formerly found experimentally by du Pont et al. [13] in the case of granular avalanche experiments in various fluids. The collapse time is shown to be similar between the dry and inertial regimes while it is longer in the viscous regime. During the early times of collapse, granular kinetic energy in the inertial and viscous regime is lower than in the dry regime. During the propagating phase, the lubrication force between the grains tends to locally decrease the intergranular friction which promotes spreading while the viscous drag slows down the grains. Overall, it turns out that the length of final deposits in the dry regime is larger than in the viscous but could be similar in inertial regimes. These numerical simulations highlighted the complex influence of the fluid phase on the collapse dynamics. It has to be noted that the influence of the initial packing fraction has not yet been reported numerically. Moreover, due to the local nature of the pore pressure feedback at the grain scale, the quantitative influence of the packing fraction should be addressed by three-dimensional simulations. Here, we study the collapse of granular columns in a fluid by three-dimensional simulations in the viscous regime. An immersed boundary method (IBM) which resolves the fluid flow around moving non-deformable particles is coupled to a soft-sphere discrete element method
(DEM) which computes the Lagrangian motion of the particles, including grain-grain interactions. A lubrication force is added to the equation of the grains motion to properly capture rebound interactions in a fluid.

This article is structured as follows. First, the description of the numerical methods is presented. Then, the three-dimensional collapse of granular columns in a fluid is simulated in the viscous regime. After some preliminary tests, we pay a special attention to the effect of the initial packing fraction and aspect ratio on the collapse dynamics.

# 2 Numerical Methods

In this section, the coupled IBM/DEM is outlined. The reader may refer to [11] for a detailed description and validation of the numerical approaches and coupling.

## 2.1 Fluid Flow Computation

The Immersed Boundary Method [14], denoted by IBM, solves the modified Navier-Stokes equations (1) and (2) for a Newtonian fluid on a Cartesian grid in the entire computational domain, e.g. in the fluid domain as well as inside the particles.

$$ \nabla \cdot \mathbf{V} = 0, \qquad (1) $$

$$ \rho \frac{\partial \mathbf{V}}{\partial t} + \rho \mathbf{V} \cdot \nabla \mathbf{V} = \rho \mathbf{g} + \nabla \cdot [\mu(\nabla \mathbf{V} + t'\nabla \mathbf{V})] - \nabla P + \rho \mathbf{f}, \qquad (2) $$

where $\rho$ and $\mu$ denote the fluid density and dynamic viscosity whereas $V$ and $P$ are the local fluid velocity and pressure, respectively, and $\mathbf{g}$ is the gravity vector. A body-force source term, namely $\mathbf{f} = \alpha(U - V)/\Delta t_f$, is added to the Navier-Stokes momentum equation (2) so that the presence and motions of the solid particles in the fluid are accounted for. $U$ and $\alpha$ are the local velocity of the solid particle and the solid volume fraction, respectively, and $\Delta t_f$ is the fluid time step. The hydrodynamic force and torque generated by the fluid upon the particle $p$ of volume $\vartheta_p$ are calculated as

$$ \mathbf{F}_h = - \frac{\rho \rho_p}{\rho_p - \rho} \int_{\vartheta_p} f \mathrm{d}V, \qquad (3) $$

$$ \tau_h = - \frac{\rho \rho_p}{\rho_p - \rho} \int_{\vartheta_p} \mathbf{r} \times f \mathrm{d}V, \qquad (4) $$
where **r** is the position vector inside the solid particle, starting from the centre of mass and $\rho_p$ is the particle density.

## 2.2 Particle Motion Computation

A discrete element method is used to compute the dynamics of the granular phase, described as a collection of spherical particles in contact. Newton's equations (5) and (6) for the particle linear and angular momentum, respectively, are solved to calculate the motion of a solid element *p*, of mass $m_p$, evolving in a collection of grains entirely immersed in a fluid:

$$m_p \frac{du_p}{dt} = m_p g + \sum_{j \neq p} F_{c-pj} + \sum_{j \neq p} F_{l-pj} + F_h, \quad (5)$$

$$I_p \frac{d\omega_p}{dt} = \sum_{j \neq p} \tau_{c-pj} + \tau_h, \quad (6)$$

where $u_p$ and $\omega_p$ are, respectively, the translational and rotational velocities of the particle *p*, and $I_p = \frac{2}{5} m_p R_p^2$ is the inertia coefficient of the spherical particle *p* of radius $R_p$. The forces applied on the particle *p* can be divided in four types of contributions, namely the particle weight; the contact forces $F_{c-pj}$ which arise from the contact between particles *p* and *j*, the lubrication forces $F_{l-pj}$ which emerge from the liquid film drainage when other particles *j* move closely to the particle *p* and the hydrodynamic force $F_h$, defined in (3). The torques applied to the particle *p* originate from contact interactions of other particle *j*, e.g. $\tau_{c-pj}$, and from the local hydrodynamics near particle *p*, e.g. $\tau_h$ defined in (4). We define these force and torque contributions in the following.

Solid-solid interactions are modelled with a DEM method [15]. We consider the normal unit vector **n** which links the centre of mass of the two solids involved in a contact interaction. A tangent unit vector **t** to **n** is defined and permits to decompose the contact force such as $F_{c-pj} = F_{c-pj}^n \mathbf{n} + F_{c-pj}^t \mathbf{t}$. A small overlap, noted $\delta_n$, between the solid spheres during the contact allows to determine the normal contribution $F_{c-pj}^n$ of the contact force with a damped mass-spring model as

$$F_{c-pj}^n = \begin{cases} 0 & \text{if } \delta_n > 0 \\ \max(0, -k_n \delta_n - \gamma_n \frac{d\delta_n}{dt}) & \text{otherwise,} \end{cases} \quad (7)$$

where $k_n$ and $\gamma_n$ are the normal stiffness and the damping coefficient, respectively, and are defined as a function of the dry restitution coefficient $\epsilon_{max}$ and the contact time $t_c$ by $\gamma_n = -2m_* \ln(\epsilon_{max})/t_c$ and $k_n = m_* \pi^2 / t_c^2 + \gamma_n^2 / 4m_*$, where $m_* = m_p m_j / (m_p + m_j)$ is the effective mass.
The tangential contribution $F_{c-pj}^n$ of the contact force is computed with a mass-spring model while satisfying the Coulomb's law as follows

$$F_{c-pj}^t = -\min(|k_t \delta_t|, |\mu_c F_{c-pj}^n|) \operatorname{sign}(\delta_t), \quad (8)$$

$$\frac{d\delta_t}{dt} = (\mathbf{u}_p - \mathbf{u}_j) \cdot \mathbf{t} \qquad (9)$$

where $\mu_c$ is the friction coefficient between the particle $p$ and $j$, $\delta_t$ is the tangential penetration which verifies at all times the differential equation (9) starting from 0 at the beginning of the contact interaction. The tangential stiffness is defined by $k_t = 0.2k_n$ as in [16].

The torque due to tangential contact interaction between the particle $p$ and $j$ is calculated as $\tau_{c-pj} = R_p \mathbf{n} \times F_{c-pj}^t \mathbf{t}$.

When two particles approach or separate from each other in a fluid, the flow structure may not be entirely captured by the present IBM method. As shown by several authors (see [17–20], for example), this issue is addressed using a normal lubrication force $F_{l-pj} = F_{L-pj} \mathbf{n}$ which reads

$$F_{L-pj} = \begin{cases} -\frac{6\pi\mu(\mathbf{u}_p \cdot \mathbf{n} - \mathbf{u}_j \cdot \mathbf{n})}{\delta_n + \eta_e} R_*^2 & \text{if } 0 \le \delta_n \le \frac{R_*}{2} \\ 0 & \text{otherwise,} \end{cases} \quad (10)$$

where $R_* = R_p R_j / (R_p + R_j)$ is the effective radius of the particles $p$ and $j$ and $\eta_e$ is a so-called effective roughness length which accounts for particle surface asperities and may range in $[10^{-6}R_*; 10^{-3}R_*]$ (see [19]). A contact time step $\Delta t_c = t_c/100$ is set to resolve the time integration of (5), (6) and (9).

## 2.3 Particles-Fluid Numerical Coupling

A numerical coupling of the immersed boundary method with the discrete element method is here done based on the assumption that the fluid timescale $t_f$ is much larger than the contact time $t_c$. Consequently, one can solve the particle motion equations (5) and (6) with (9) separately from the Navier–Stokes equations (1) and (2) and considering that the fluid is frozen on a time scale of the order of $t_f$.

The fluid pressure and velocity are first calculated with the IBM in the whole domain, second, the contributions of the hydrodynamic forces and torques on each particles in the system are computed and transmitted to the DEM solver. Then, the new position and velocity of particles are computed taking into account contact and hydrodynamics interactions including lubrication, using several contact time steps within a fluid time step. Subsequently, the new particles position and velocity are given back to the IBM solver.
### 3 Collapse of a Granular Column in a Fluid

The three-dimensional IBM/DEM simulations of the collapse of a granular column, composed of approximately $O(10^3)$ particles, are presented for a Stokes number $St = \frac{\sqrt{\rho_p D \Delta \rho g D}}{18\sqrt{2}\mu} \approx 1.7$, a density ratio $r = \sqrt{\frac{\rho_p}{\rho}} = 2.8$, a Reynolds number $Re = \frac{St}{r} \approx 0.6$ and a Galileo number $Ga = \rho \sqrt{(\frac{\rho_p}{\rho} - 1)gD^3}/\mu \approx 15.5$. In the view of the avalanche regime classification of Courrech [13] in the $(St, r)$ parameter space, one can expect our simulations to be in the viscous regime.

Throughout this work, the gravity is set to $9.81 \text{ m/s}^2$ and is pointing in the $-y$ direction (see the Cartesian system of reference in Fig. 1), solid contact parameters are $\epsilon_{max} = 0.87$, $t_c \sqrt{g/D} = 2.5 \times 10^{-3}$ and $\mu_c = 0.25$, the effective roughness length in the lubrication model is set to $\eta_c/R_* = 2.10^{-4}$, and the fluid parameters are $\mu = 1 \text{ Pa s}$, $\rho = 1000 \text{ kg/m}^3$ and $\Delta t_f \sqrt{g/D} = 2.5 \times 10^{-2}$. To avoid crystallization phenomena in the granular medium, the diameter of particles is varied in a range of $\pm 5\%$ from its mean $D = 15.2 \text{ mm}$. A no-slip boundary condition is imposed on all walls of the domain except in the spanwise direction $z$,

<figure><img src="image_2.png" /><figcaption>**Fig. 1** Initial configurations used for simulations. **a** 1 and 2, **b** 3, **c** 4 and **d** 5. The Cartesian system of reference $(x, y, z)$, the size of the fluid domain $[L: H; W]$ and the initial dimensions of the granular column $[R_o; H_o]$ are defined. Half spheres are glued to the bottom wall to imitate a rough base</figcaption></figure>
where periodic boundary conditions are applied. Also, half-grains are glued on the bottom wall in order to numerically mimic a rigid rough bottom as seen in Fig. 1. The solid particles of the bottom roughness have the same contact parameters than the ones free to move. All the simulations performed are summarized in Table 1 with the corresponding geometrical and numerical parameters.

In the following, some preliminary observations and tests on the fluid grid resolution $\Delta x$ and on the spanwise length $W$ of the column are discussed. Then, a simulation of collapse is presented in details from its initial configuration to its final state where all gains are deposited and statics. After, we investigate the effect of the initial aspect ratio and packing fraction on the granular dynamics and deposits.

# 3.1 Preliminary Observations

## 3.1.1 Discussion About the Initial Packing Fraction

The initial granular columns can be obtained using only the DEM approach by letting the particles freely fall in an elongated box which permits to create granular columns with an important initial packing fraction $\phi_i = 0.6 \pm 0.002$. This packing fraction is measured inside the granular column for simulations 1–5 where the solid fraction $\alpha$, issued from the IBM, is averaged in a given control volume. A 40% variation of the control volume leads to a 0.7% variation of the packing fraction. As done in granular collapse experiments of [6], another way to estimate $\phi_i$ is to calculate the ratio between the volume occupied by the grains $1/\rho_p \sum_{i=1}^{N_g} m_i$ and the granular column volume $R_o H_o W_o$. In our case, this estimation leads to a packing fraction 0.585 which is 2.5% less than the calculation with the solid fraction of the IBM. This subtle but noticeable difference may be due to the definition and estimation of the average height of the column $H_o$. We shall come back to this point in the results section. Note that, this initial state will be referred as the dense packing state in the following. In order to make a granular column with a loose state, the viscous dissipation on the particle dynamics during sedimentation is accounting for by using the IBM/DEM coupling. In that case, the viscous drag modifies the dynamics of the grains and allows an initial packing fraction $\phi_i = 0.58$ that we will denote as the loose case (see simulation 6 in Table 1). A similar 2.6% bias on this estimation is obtained when measuring the ratio of grains volume and box volume, as in the dense case.

## 3.1.2 A Typical Granular Collapse

IBM/DEM simulation of a granular collapse for $a = 0.5$, $\phi_i = 0.6$ and $D/\Delta x = 10$, is considered here with a spanwise direction size $5D$ (see simulation 3 in Table 1). The mesh size used for simulation 3 is $[N_x; N_y; N_z] = [512; 176; 52]$ and the initial column is made of 798 grains. The meshes are uniform in size where the granular
<table><caption>**Table 1** Parameters used is the present simulations: the initial aspect ratio $a = H_o/R_o$, and packing fraction $\phi_i$, the number of grains $N_g$, the cells number in each direction in the domain $[N_x; N_x; N_z]$ as well as the size of the fluid domain $[L; H; W]$ presented in Fig. 1, the fluid grid resolution $\Delta x$, the initial dimensions of the granular column $[R_o; H_o; W_o]$ and the maximum horizontal kinetic energy averaged over all grains, noted $E_{kx}^{\text{max}}$, during the simulation are presented</caption><thead><tr><th>Simulation</th><th>a</th><th>φ<sub>i</sub></th><th>N<sub>g</sub></th><th>N<sub>x</sub>; N<sub>y</sub>; N<sub>z</sub></th><th>L; H; W</th><th>Δx</th><th>R<sub>o</sub>; H<sub>o</sub>; W<sub>o</sub></th><th>E<sub>kx</sub><sup>max</sup>/(mgD)</th></tr></thead><tbody><tr><td>1</td><td>0.5</td><td>0.6</td><td>327</td><td>512; 176; 22</td><td>66D; 30D; 2D</td><td>D/10</td><td>16D; 8D; 2D</td><td>0.023</td></tr><tr><td>2</td><td>0.5</td><td>0.6</td><td>327</td><td>1008; 344; 42</td><td>66D; 30D; 2D</td><td>D/20</td><td>16D; 8D; 2D</td><td>0.020</td></tr><tr><td>3</td><td>0.5</td><td>0.6</td><td>798</td><td>512; 176; 52</td><td>66D; 30D; 5D</td><td>D/10</td><td>16D; 8D; 2D</td><td>0.015</td></tr><tr><td>4</td><td>1</td><td>0.6</td><td>665</td><td>800; 288; 27</td><td>98D; 46D; 2D</td><td>D/10</td><td>16D; 16D; 2D</td><td>0.082</td></tr><tr><td>5</td><td>1.5</td><td>0.6</td><td>985</td><td>688; 352; 22</td><td>88D; 54D; 2D</td><td>D/10</td><td>16D; 24D; 2D</td><td>0.452</td></tr><tr><td>6</td><td>0.5</td><td>0.58</td><td>318</td><td>512; 176; 22</td><td>66D; 30D; 2D</td><td>D/10</td><td>16D; 8D; 2D</td><td>0.032</td></tr></tbody></table>
<figure><img src="image_3.png" /><figcaption>**Fig. 2** Three-dimensional IBM/DEM simulation of a granular column collapse composed of approximately 800 grains in a viscous fluid. Iso-contours of pressure are represented in a vertical plane perpendicular to the spanwise direction. The fluid pressure scaled by $\Delta\rho gH_o$ (where $\Delta\rho = \rho - \rho_p$) between consecutive iso-contours is $9 \times 10^{-3}$. An iso-surface of vorticity is plotted for a value of 0.4 and is represented in light blue. At $t = 0^{-}$, a wall which supports the granular column is plotted in dark blue plane. The grains are represented by an iso-surface of solid fraction $\alpha = 0.8$. Time and vorticity are scaled by $\sqrt{D/g}$, $\sqrt{g/D}$, respectively</figcaption></figure>

media may flow, e.g. in the region $0 \le x/D \le 50$, $0 \le y/D \le 16$ and $0 \le z/D \le 5$; whereas the grid size outside of this region is varied following an arithmetic progression up to walls. This permits to put the limits of the domain far enough so that it does not influence the collapse dynamics and it saves computational resources. Initially, a wall represented in dark blue in Fig. 2 at the time $t\sqrt{g/D} = 0^{-}$ and supporting the granular column, is released so that the collapse occurs in the $x$ direction. Several instants during the collapse are depicted in Fig. 2 where iso-surface of vorticity is plotted as well as iso-contours of pressure in a vertical plane perpendicular to the spanwise direction. Initially, the pressure is in hydrostatic equilibrium and the fluid and grains are at rest. At the time $t\sqrt{g/D} = 12.7$, vorticity stemming from a local shear stress is created in the avalanche zone, e.g. where the grains flow. At this moment, the avalanche front velocity is maximum as seen in Fig. 4b. Then, for $25.4 \le t\sqrt{g/D} \le 38.1$, the granular front abruptly slows down and residual vorticity at the surface of the granular deposit is dissipated while the pressure relaxes towards the hydrostatic equilibrium. In the present simulation, the granular deposit has a final length $25D$ and a final high $8D$ forming a trapezoidal form.
### 3.1.3 Effect of the Spatial Resolution

Two simulations with $a = 0.5$, $\phi_1 = 0.6$, $W = 2D$ and $N_p = 327$ are considered for two grid resolutions, namely $\Delta x = D/10$ in simulation 1 and $\Delta x = D/20$ in simulation 2. The meshing used for simulation 1 (resp. 2) is $[N_x; N_y; N_z] = [512; 176; 22]$ (resp. $[1008; 344; 42]$). The time evolution of the outer envelope of the granular column is plotted for simulations 1 and 2 in Fig. 3. Also, the granular front location is plotted as a function of time for these simulations in Fig. 4a. The temporal evolution of the front is similar up to $t\sqrt{g/D} = 17$. Then, its position seems frozen for simulation 2 while it continues to move for simulation 1, i.e. the dissipation due to the fluid is slightly underestimated for a lower spatial resolution of the fluid phase. When the grains are statics, the granular deposit has a similar shape in both cases. The relative position of the front between simulation 1 and 2 is 10%. Although this 10% error in the final front position is somehow not negligible, the time evolution of the front position highlights a similar trend in both cases. Given the computational cost of the highest spatial resolution, a spatial resolution

<figure><img src="image_4.png" /><figcaption>**Fig. 3** Temporal evolution of the outer envelope of the granular column during the collapse in a fluid with a spatial resolution. **a** $\Delta x = D/10$ and **b** $\Delta x = D/20$. The time is scaled by $\sqrt{D/g}$ and varies from 0 to 37.2 with steps of 2.48. Note that the collapse starts at time $t = 0$</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>**Fig. 4** Temporal evolution of the granular front location scaled by the initial length of the column (($x_{\text{front}} - R_o)/R_o$ **a** o (resp. □) for simulation 1 (resp. 2) with a spatial resolution $\Delta x = D/10$ (resp. $\Delta x = D/20$) and **b** o (resp. x) for simulation 1 (resp. 3) with a spanwize direction size $W = 2D$ (resp. $W = 5D$)</figcaption></figure>
$\Delta x = D/10$ seems adapted to perform simulations of granular collapse in the configuration studied here. In the following, we will, therefore, use this spatial resolution.

### 3.1.4 Effect of the Domain Size in the Spanwise Direction

Two simulations with $a = 0.5$ and $\phi_i = 0.6$ are considered for two different spanwise directions $W = 2D$ (simulation 1) and $W = 5D$ (simulation 3). The granular front is plotted as a function of time in Fig. 4b. In both cases, the spreading is similar and the granular deposit has the same shape (not shown here). Nevertheless, we note that the collapse in simulation 3 is faster than in simulation 1 around $t\sqrt{g/D} = 7$. In addition, at $t\sqrt{g/D} = 20$, the granular front is static in simulation 1 while it continues to slowly move in simulation 3 up to $t\sqrt{g/D} = 40$ where the static state is reached. These differences may be a consequence of the granular confinement in the spanwise direction which tends to modify the rearrangement of grains at low speed as it is the case in the present simulation in the viscous regime and for a small aspect ratio. Note that, the quantitative influence of the confinement observed here is of the same order as the one due to the spatial resolution discussed in the previous section.

The uncertainties on the granular dynamics should, therefore, be accounted for but are actually small compared with the influence of the aspect ratio in the range considered in this study. In the light of these preliminary tests, we, therefore, decide to perform a parametric study of the granular collapse in a fluid with a spatial discretization $\Delta x = D/10$ and a computational domain width of $2D$.

## 3.2 Effect of Initial Aspect Ratio

Three granular collapses are simulated with the IBM/DEM approach with aspect ratios $a = [0.5; 1; 1.5]$ and an initial packing fraction $\phi_i = 0.6$ (simulations 1, 4 and 5 in Table 1). We recall that the fluid and particles properties are the same than presented in previous sections.

Following Courrech du Pont et al. [13], the time scale of the granular collapse in a fluid could be controlled by various temporal scaling, namely the characteristic viscous time $T_v = \frac{18\mu}{\Delta\rho g H_o}$, inertial time $T_i = \sqrt{\frac{\rho H_o}{2\Delta\rho g}}$ and dry or free-fall time $T_{ff} = \sqrt{\frac{2\rho_p H_o}{\Delta\rho g}}$. As already explained at the beginning of Sect. 3, in the view of the avalanche regime classification of Courrech [13] in the $(\text{St}, r)$ parameter space, our collapse simulations are expected to be in the viscous regime, meaning that the motion of grains is controlled by the viscous time $T_v$. For the presented simulations, the time evolution of the front is plotted in Fig. 5, where the time is scaled by $T_v$, $T_i$ and $T_{ff}$. The front position accelerates first then the front velocity is approximately
<figure><img src="image_6.png" /><figcaption>**Fig. 5** Position of the granular front ($x_{front} - R_o/R_o$) with different initial aspect ratios $a = [0.5; 1; 1.5] = [\Delta; \bigcirc; \Box]$ versus time scaled by $T_v$, $T_i$ and $T_{ff}$. In all cases, the front propagation has an S-form where three phases are observed: the acceleration, the constant velocity propagation and the deceleration. A red dashed line in the left graph shows a master curve for a large portion of the propagation phase for all simulations, which shows the relevance of the viscous time scale $T_v$ in the present simulations for all aspect ratios</figcaption></figure>

constant and, at last, it decelerates abruptly for all studied aspect ratios. When the viscous time $T_v$ is used to scale the time, a linear master curve is found in Fig. 5 (on the left) which corresponds to a constant velocity for the front propagation equals to approximately $10^{-3}R_o/T_v$ m/s. This confirms that the dynamics of the dense granular flow generated by the collapse corresponds to a viscous regime in the range of parameters considered here.

Let us consider, the total granular horizontal kinetic energy defined as $E_{kx} = \sum_{i=1}^{N_g} \frac{1}{2} m_i u_{ix}(t)^2$ where $u_{ix}$ is the streamwise velocity of grain $i$ and $m_i$ is the mass of grain $i$. The maximal horizontal kinetic energy for the whole duration of a simulation is noted $E_{kx}^{\max}$ and is given in Table 1. It can be noted that $E_{kx}^{\max}$ increases with the initial aspect ratio $a$ as this energy is transferred at early stages of the collapse from the initial potential energy of the column which increases with $a$. In the case $a=1$, Topin et al. [12] found in their two-dimensional numerical simulations $E_{kx}^{\max}/(mgD) \approx 0.02$ whereas we get 0.082 for the same $a$ (simulation 4). Even if these values are of the same order of magnitude, a factor 4 exists. This quantitative discrepancy between 2D and 3D simulations would deserve a deeper investigation. In the current state of research, one can conclude that 3D simulations could enhance the horizontal kinetic energy of the granular mass compared with 2D simulations.

The final granular deposit height $H_\infty$ and length $R_\infty$ are compared to the experiments of [6] and plotted in Fig. 6. The final deposit length obtained from the IBM/DEM simulations is in accordance with the experimental measurements for a similar initial packing fraction. One can note that the final deposit height found in the simulations for $a \ge 1$ lies in between the experimental results of Rondon [6], where the initial packing fraction varies from 0.55 to 0.6. Recalling the discussions
<figure><img src="image_7.png" /><figcaption>**Fig. 6** Final deposit height $H_\infty/R_o$ (left) and final deposit length $R_\infty/R_o$ (right) as a function of the aspect ratio $a$. □: IBM/DEM simulations for $\phi_i = 0.6$; ○ and ●: experiments of Rondon [6] for an initial packing fraction $\phi_i = 0.55$ and $\phi_i = 0.6$, respectively; ◇: 2D simulations of Topin et al. [12] in the viscous regime for an initial packing fraction $\phi_i = 0.8$</figcaption></figure>

in Sect. 3.1.1 about the discrepancy estimation on the initial packing fraction, this numerical prediction seems consistent with experiments. For comparison, two-dimensional numerical predictions of Topin et al. [12] are represented in Fig. 6. As a conclusion, in the viscous regime, it seems that there is a quantitative agreement between the 2D numerical simulations of [12] and the 3D present simulations regarding the morphology of the deposits.

## 3.3 Effect of Packing Fraction

Rondon et al. [6] have observed in their experiments an effect of the initial packing fraction $\phi_i$ on the granular collapse dynamics in the viscous regime. At the end of the collapse, all grains are static and, as shown by [6], the granular deposit can be found in a triangular or trapezoidal form depending on the initial aspect ratio $a$ and the initial packing fraction $\phi_i$. The geometry of the granular deposit for the experiments of [6] and the present IBM/DEM simulations are depicted in the $(\phi_i, a)$ parameter space in Fig. 7. A good qualitative agreement is found between experiments and present simulations. In particular, the shape evolution from trapezoidal to triangular shape for decreasing $\phi_i$ at small $a$ and increasing $a$ observed in the experiments is also obtained in the simulations. However, a quantitative discrepancy is obtained as the simulation 6, i.e. for $a = 0.5$ and $\phi_i = 0.58$, leads to a triangular form whereas [6] found a trapezoidal one for this specific value of $(\phi_i, a)$. This discrepancy on the transition from triangular to trapezoidal shape as a function of $\phi_i$ may be explained by the procedure to estimate the initial packing fraction $\phi_i$ in the simulation which is different from the one used in the experiments [6], as mentioned in Sect. 3.1.1. These definitions of the packing fraction may not be equivalent and an artificial offset is likely to exist between those values. For example, a 4% offset in the direction of the small values of $\phi_i$ would make all the present numerical results to be in full agreement with the experiments of [6]. Recall that computing the initial packing fraction from the local integration of the solid
**Fig. 7** Granular deposit forms observed in the presented simulations (△ for triangle and □ for trapeze) and in the experiments of Rondon et al. ([6]) (▲ for triangle and ■ for trapeze) in the ($\phi_i, a$) parameter space

<figure><img src="image_8.png" /></figure>

**Fig. 8** Time evolution of granular fronts for the dense (○) and the loose (♦) cases; simulations 1 and 6 respectively. Time is scaled by the viscous time $T_v$ (see the definition in the text in Sect. 3.3)

<figure><img src="image_9.png" /></figure>

volume fraction or from the average height of the granular leads to a 2.5% difference from our simulation results.

Let us define, a reference case which is simulation 1 with $a = 0.5$ and $\phi_i = 0.6$ (see Table 1) that we will refer to the dense case. The temporal evolution of the granular front for the loose ($\phi_i = 0.58$) and the dense ($\phi_i = 0.6$) cases are plotted in Fig. 8. The collapse in the loose case occurs earlier than in the dense case with a delay of approximately 1000$T_v$. Also, the final deposit length in the loose configuration is greater by 10% than the dense one. This has been already found in the experimental work of [6] and confirms the influence of the initial packing fraction on the collapse dynamics of a granular column immersed in a viscous fluid.

Several instances of these two collapses (loose and dense) are depicted in Fig. 9 showing the pressure field which has been averaged in the spanwise direction, and the outer envelope of the granular column. At early times, e.g. at $t = 900$ (scaled by the viscous time $T_v$), a low-pressure (resp. high-pressure) region, relative to the hydrostatic pressure, develops inside the column in the dense (resp. loose) configuration. This is a consequence of the rearrangement of the granular media which
<figure><img src="image_10.png" /><figcaption>**Fig. 9** Time evolution of the fluid pressure with the form of the external envelop for two collapses of aspect ratio of *a* = 0.5, Δ*x* = *D*/10 and *φ*ᵢ = 0.6 on the left and *φ*ᵢ = 0.58 on the right. The time is scaled by *T*_v, defined in Sect. 3.2. The pressure fields are averaged over the spanwise direction</figcaption></figure>
undergoes a decompaction and compaction, respectively. At the time $t = 2700$, the low-pressure region is always present in the dense case while it appears in the loose case as a consequence of the granular column motion. At the time $t = 3600$, granular fronts of the two columns are spreading and the low-pressure zones tend to disappear while a high-pressure is observed at the front of the granular avalanche in both cases. Once the column is at rest, pressure relaxes up to the hydrostatic equilibrium. At time $t = 18000$, the final deposit has a triangular (resp. trapezoidal) form in the loose (resp. dense) case. The measured static avalanche angles are approximately $26^\circ$ in both cases. These results on the final form of the granular deposits are in qualitative accordance with the experiments of Rondon et al. [6].

The initial packing fraction has an influence on grains rearrangement during the granular collapse in a viscous fluid. In the loose state, the fluid initially located in the pores is pushed away from the column by the grains which leads to the creation of high-pressure regions. In this case, the fluid flow from the inner region of the granular medium to the outer region can occur simultaneously with the granular collapse. On the other hand, for an initially dense granular column, a local decompaction of the granular medium is necessary prior any motion of the grains, leading to low-pressure regions. This low-pressure region is associated with a fluid flow from the outer region into the inner region. The fluid flow occurs on a viscous time scale and delays the granular collapse phenomenon. To quantify this pressure phenomena and to compare to the experiments, the time evolution of the basal pore pressure, noted $P_b$ and defined as the pressure under the column (at a distance $x = 0.36R_0$ and $x = 0.8R_0$ in the experiments and simulations, respectively), is plotted in Fig. 10 for the experiments of Rondon et al. [6] (a) and present simulations (b). We observe a similar trend for the two methods: an initial zero value for basal pressure is followed by a negative or positive peak pressure depending on the initial packing fraction, and eventually returns to zero relative pressure at time $t = 15000T_v$. Note that the pressure signals between the dense and the loose cases look symmetrical in both the simulations and the experiments. Nevertheless, the dispersion of the pressure signal obtained by numerical simulations is greater than in the experiments, which is to be expected since the number of grains involved in the collapse is smaller in the simulations. One may define the signed maximal pressure over time, namely $P_b^{\max}$ which is plotted versus the initial packing fraction $\phi_i$ in Fig. 10c for the experiments of Rondon et al. [6] and the present simulations. Quantitative agreement is observed, which confirms that the pore pressure feedback phenomenon observed in a viscous regime is captured by the IBM/DEM coupling model.

# 4 Conclusions

The three-dimensional unsteady collapse of the granular column in a viscous fluid has been investigated with an IBM/DEM approach. Present numerical simulations allow one to confirm quantitatively several experimental observations of Rondon et al. [6] regarding morphology, characteristic sizes of granular deposits and the
<figure><img src="image_11.png" /><figcaption>**Fig. 10** a Time evolution of the basal pore pressure $P_b$ (relative to the hydrostatic pressure) in the granular column collapse in some fluid experiments of Rondon et al. [6] for an initial aspect ratio $a = 0.67$ for a dense (●) and a loose (○) initial configurations where the initial packing fractions are $\phi_i = 0.6$ and $\phi_i = 0.55$, respectively. This image is extracted from [6] and the probe is located at $x = 0.36R_0$. b Temporal evolution of the basal pressure (at a distance $x = 0.8R_0$) scaled by $\Delta\rho gH_o$ in the present IBM/DEM simulations for a dense (●) and a loose (○) initial cases with $\phi_i = 0.6$ and $\phi_i = 0.58$, respectively. c Maximal basal pressure $P_b^{\max}$ versus initial packing fraction $\phi_i$: experiments of Rondon [6] are represented with small black squares (■) and present IBM/DEM simulations for a dense (●) and a loose (○) initial granular states</figcaption></figure>

basal pressure below the column. In the presented simulations, the collapse dynamics is controlled by the viscous time $T_v$. To our knowledge, a numerical approach, e.g. the IBM/DEM method, is able for the first time to capture the pore pressure feedback phenomenon in flowing fluid-grains mixture. The effect of the initial packing fraction has a great influence of the dynamics of granular collapse in the simulation results as in the experiments of Rondon et al. [6]. Furthermore, the IBM/DEM permits to investigate the inner state of the granular column during the collapse, in particular, the evolution of the pressure field inside the granular column can be analyzed which is difficult to do in experiments. Simulations of the collapse of a granular column immersed in a fluid can be performed in the inertial and free-fall regimes as well varying the nature of the fluid and/or the particles. This opens the way to perform a parametric study on the collapse dynamics in the parameter spaces ($\rho_p/\rho, St$) and ($\phi_i, a$). Furthermore, local information from these simulations may give a better understanding of the rheology of these complex grains-fluid mixture dynamics; these numerical perspectives should be investigated in a future work.

**Acknowledgements** This work was granted access to the HPC resources of CALMIP supercomputing center under the allocation P1027, the contributions of which is greatly appreciated. The study has been supported by the ‘Agence Nationale de la Recherche’ in the frame of the project ANR-12-2013-ModSed.
References

1. Lube, G., Huppert, H. E., Sparks, R. S. J., & Hallworth, M. A. (2004). Axisymmetric collapses of granular columns. *Journal of Fluid Mechanics*, **508**, 175–199.

2. Lajeunesse, E., Mangeney-Castelnau, A., & Vilotte, J. P. (2004). Spreading of a granular mass on a horizontal plane. *Physics of Fluids*, **16**(7), 2371–2381 (1994-present).

3. Lube, G., Huppert, H. E., Sparks, R. S. J., & Freundt, A. (2005). Collapses of two-dimensional granular columns. *Physical Review E*, **72**(4), 041301.

4. Balmforth, N. J., & Kerswell, R. R. (2005). Granular collapse in two dimensions. *Journal of Fluid Mechanics*, **538**(1), 399–428.

5. Thompson, E. L., & Huppert, H. E. (2007). Granular column collapses: Further experimental results. *Journal of Fluid Mechanics*, **575**, 177–186.

6. Rondon, L., Pouliquen, O., & Aussillous, P. (2011). Granular collapse in a fluid: role of the initial volume fraction. *Physics of Fluids*, **23**(7), 073301.

7. Iverson, R. M. (2000). Landslide triggering by rain infiltration. *Water Resources Research*, **36** (7), 1897–1910.

8. Zenit, R. (2005). Computer simulations of the collapse of a granular column. *Physics of Fluids*, **17**(3), 031703 (1994-present).

9. Staron, L., & Hinch, E. J. (2005). Study of the collapse of granular columns using two-dimensional discrete-grain simulation. *Journal of Fluid Mechanics*, **545**, 1–27.

10. Girolami, L., Hergault, V., Vinay, G., & Wachs, A. (2012). A three-dimensional discrete-grain model for the simulation of dam-break rectangular collapses: Comparison between numerical results and experiments. *Granular Matter*, **14**(3), 381–392.

11. Izard, E., Bonometti, T., & Lacaze, L. (2014). Simulation of an avalanche in a fluid with a soft-sphere/immersed boundary method including a lubrication force. *The Journal of Computational Multiphase Flows*, **6**(4), 391–405.

12. Topin, V., Monerie, Y., Perales, F., & Radjai, F. (2012). Collapse dynamics and runout of dense granular materials in a fluid. *Physical Review Letters*, **109**(18), 188001.

13. du Pont, S. C., Gondret, P., Perrin, B., & Rabaud, M. (2003). Granular avalanches in fluids. *Physical Review Letters*, **90**(4), 044301.

14. Yuki, Y., Takeuchi, S., & Kajishima, T. (2007). Efficient immersed boundary method for strong interaction problem of arbitrary shape object with the self-induced flow. *Journal of Fluid Science and Technology*, **2**(1), 1–11.

15. Cundall, P. A., & Strack, O. D. (1979). A discrete numerical model for granular assemblies. *Geotechnique*, **29**(1), 47–65.

16. Shäfer, J., Dippel, S., & Wolf, D. E. (1996). Force schemes in simulations of granular materials. *Journal de Physique I*, **6**(1), 5–20.

17. Kempe, T., & Fröhlich, J. (2012). Collision modelling for the interface-resolved simulation of spherical particles in viscous fluids. *Journal of Fluid Mechanics*, **709**, 445–489.

18. de Motta, J. B., Breugem, W. P., Gazanion, B., Estivalezes, J. L., Vincent, S., & Climent, E. (2013). Numerical modelling of finite-size particle collisions in a viscous fluid. *Physics of Fluids*, **25**(8), 083302 (1994-present).

19. Izard, E., Bonometti, T., & Lacaze, L. (2014). Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid. *Journal of Fluid Mechanics*, **747**, 422–446.

20. Costa, P., Boersma, B. J., Westerweel, J., & Breugem, W. P. (2015). Collision model for fully resolved simulations of flows laden with finite-size particles. *Physical Review E*, **92**(5), 053012.