E3S Web of Conferences 415, 02025 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502025

Micromechanical simulations of the collapse of a submerged
granular column
Francesco Zarattini1*, Antonio Pol2, Fabio Gabrieli1 and Bruno Chareyre3
1University of Padova, ICEA dept., Italy
2Université Gustave Eiffel, MAST/GPEM, France
3Université Grenoble Alpes, 3SR, G-INP, CNRS, FranceAbstract. In this work, we present a numerical study of a granular collapse in a viscous fluid adopting aDEM-PFV (Pore Finite Volume) coupled model. We analyze the time evolution of the collapse and the final
morphology of the deposit (triangular or trapezoidal) as a function of the initial solid fraction and column
aspect ratio. Finally, we discuss the role of the initial packing fraction of the granular material on the time
evolution of the pore pressure. The numerical results are consistent with experimental observation [3]showing the capability of the DEM-PFV coupled model to simulate saturated granular media from triggering
to flowing and stopping conditions.

1 Introduction
Wet granular media are frequently encountered ingeophysical problems (e.g. debris flows, landslides). As
known from classical geomechanics the interaction
between the solid and fluid phases plays a fundamentalrole in the mechanical behaviour of the medium which
is mainly controlled by the saturation degree [1] and thefine content [2]. For fully saturated soils, it entails astress exchange between the phases, which is at the baseof phenomena like consolidation and liquefaction.
In the simple case of a granular collapse,
experimental studies [3, 4] have shown that even small
differences in the initial packing fraction can lead toabrupt variations of the dynamics of the collapse in thepresence of interstitial fluid. Nevertheless, information
at the microscale can hardly be obtained fromexperimental tests. In this perspective, discrete
numerical models represent a fundamental tool andseveral coupled approaches have been proposed in thetechnical literature (e.g. DEM-CFD, DEM-LBM, DEMSPH, DEM-PFV). [5]In this work, we used a discrete element method(DEM) coupled with a pore network technique (PFV) inorder to describe the saturated granular material. We
simulate the collapse of an immersed granular column(Figure 1), and we consider two different solid fractions(loose and dense) of the column as initial conditions.
The goal of the present study is to verify the capabilityof the DEM-PFV approach to reproduce the mainfeatures of the considered physical process: (i) thecollapse dynamics, (ii) the time-evolution of the porepressure. From this perspective, we compared the*numerical results with the experimental data of Rondon
et al. (2011).

Fig. 1. Test configuration: (a) view of the initial column and(b) scheme of the virtual pressure sensors position at the baseof the column.

2 DEM-PFV coupled model
We briefly describe the adopted numerical frameworkin what follows. The interested reader can refer to [6]and [7] for further details.
The solid phase is modelled through a soft-sphereapproach [8] and we assume having a viscous flow of a
Newtonian fluid within the pore space. In order toaccount for the interaction between the two phases, i.e.
poromechanical coupling, we use the DEM-PFV
numerical scheme implemented in the open-source code
YADE [9]. It should be noted that in this work, weassume the two phases to be incompressible and we do
not consider lubrication forces. The domain isdiscretized in tetrahedral regions by using a regulartriangulation, in which the vertices of the tetrahedra are
the soil particles centres; we refer to the pore spacewithin a tetrahedron as a pore. Considering two adjacentpores, their exchanged flux is linearly related to the local
pressure gradient as a consequence of the Stokesian flow

Corresponding author: francesco.zarattini@phd.unipd.it© The Authors, published by EDP Sciences. This is an open access article distributed under the terms of the Creative Commons Attribution License 4.0(https://creativecommons.org/licenses/by/4.0/).

E3S Web of Conferences 415, 02025 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502025

regime assumption. In addition, the volume balance
equation, for a pore, permits us to link the fluxes to thevariation of the pore volume. Finally, by expressing thepore-volume change in terms of particles displacements
and writing the relation between flux and pressure local
gradient, for each pore, one can obtain a linear systemof equations which, in a matrix form, reads:
KP=E𝑥˙ +Bq+Bp

collapse and not to reproduce in a quantitative manner
the experimental results of [3]. We will refer todimensionless quantities in the following analysis. As alength scale, we use the initial height of the column H0;
the time is rescaled on a viscous time scale
Tv = 18μ/ΔρgH0, while the pressure on the quantityΔρgH0.

(1)

3.1 Time evolution of the collapse

where K is the conductivity matrix, P is the pressurevector (unknown), E relates the rate of pore-volume
change to the solid particles' velocities 𝑥˙ , Bq and Bp
represent the boundary conditions in terms of imposed
flux and pressure, respectively.
At each time step, the resolution of (1) gives thediscrete field of pore pressure. From this, drag forces canbe calculated (linearly dependent on P), and these are
introduced into Newton’s second law and applied to thesolid particles.

We analyse the influence of the initial solid fraction on
the collapse dynamics at first. In this perspective, wefocus on the time evolution of the lateral profile and theposition of the centre of mass of the solid phase. The
initial samples have a reference aspect ratio AR = 2.
We display the evolution of the profile of the
granular column in Figure 2, where the profile of thecollapsed column is reported every 0.5 s. Starting froma loose sample (ϕi = 0.56) we observe a very rapid
collapse dynamics, in which a sudden reduction of theinitial height is followed by an almost instantaneous
spreading of the column. The front travels very fast andit reaches its final position after a brief time.

3 Numerical results
In this work, we mainly focus on the role of the initialsolid fraction in the dynamics of the collapse of an
immersed granular material. Therefore, we consider two
almost identical granular prismatic samples (~2500
spherical particles) that only differ by their initial solid
fraction: a dense sample (ϕi = 0.62) and a loose sample(ϕi = 0.56). The initial geometry of the column is
described by its height H0 and aspect ratio AR = H0/L0.
The numerical parameters are reported in Table 1.
We verified the validity of the hypothesis of Stokesianflow conditions by using the criterion proposed in [10].
Gravity g = -9.81m/s2 is applied along the y-direction.

(a)(b)

Table 1. Numerical model parameters.
Mean particle diameter, d [m]0.034

Solid density, ρs [kg/m3]

2600

Fluid density, ρw [kg/m3]

1000

Viscosity, μ [Pa∙s]

20

Normal contact stiffness, kn [Pa∙m]

1.72e4

Tangential contact stiffness, kt [Pa∙m]

5.17e3Fig. 2. Time evolution of the profile of the collapsing column(AR = 2) starting from an initially (a) dense sample and a (b)loose sample.

Conversely, when we start from a dense state of thegranular column (ϕi = 0.62), we observe a much slower
spreading of the column with a progressive mobilizationof the granular material. The collapse phenomenon
starts at the top of the column with the triggering of local
failures, then it progressively evolves mobilizing thelower layers. We note that the front moves much slower
than in the loose case (compare top and bottom profilesin Figure 2) and that the final runout distance issignificantly shorter in the dense case. The timeevolution of the centre of mass of the collapsing column,
which is shown in Figure 3, leads to similar
considerations: faster collapse dynamics and a larger
spreading of the granular mass when starting from aloose initial sample.

In the analysis presented here, we mainly consider
the following aspects: (i) the time evolution of theprofile, (ii) the morphology of the final deposit, (iii) thetime evolution of the pore pressure at the base of thesample. This is analogous to the experimental work ofRondon et al. [3]. It should be underlined that theobjective of this study is to show the capability of themodel to simulate the main aspects involved in the2E3S Web of Conferences 415, 02025 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502025

a loose and a dense initial configuration. In theexperimental study of Rondon et al. [3], the basal porepressure was observed to be strongly correlated to theinitial packing fraction of the column.
With this coupled particle-based approach we areable to record the evolution of the excess of porepressure during the entire collapse and evaluate its
distribution in the deformable mass. In Figure 5 we cancompare its evolution in the dense and loose sample forthe reference case with AR = 2.
(a)

(b)

Fig. 3. Time evolution of the centre of mass position zM.

3.2 Morphology of the final deposit
Experimental observations show that the final
morphology of the collapsed column may differ whenstarting from a different packing fraction [3]. In fact,
when starting from a relatively dense state both atriangular and a trapezoidal deposit can be observed
depending on the aspect ratio of the column. Large
aspect ratios lead to a triangular morphology, while forlow aspect ratios a trapezoidal deposit is observed.
Conversely, when starting from a loose sample thetriangular morphology is the dominant one (trapezoidaldeposits are observed only for very low ARs). In orderto show the capability of the DEM-PFV approach toreproduce this behaviour, we display in Figure 4 thefinal deposit obtained from the collapse of two columns
with an aspect ratio of 0.5 and 2 respectively. For bothgeometries, the initially loose (ϕi = 0.56) and dense(ϕi = 0.62) cases were considered.
(a)

(b)

Fig. 5. Time evolution of the excess of pore pressure contour(normalized) for the (a) dense and (b) loose state with AR = 2.
Red values for positive pressure (higher than the hydrostatic)
and blue areas for the negative.

We observe that a negative pore pressure develops
during the collapse phase when starting from a dense
sample. In fact, the granular material has to (locally)
dilate for the grain motion to be permitted and thisdilation causes a small suction (pore pressure lower thanthe hydrostatic one) of the fluid phases from the outerregion (i.e., p < 0). Accordingly to Terzaghi’s stressprinciple, the development of a negative excessive porepressure translates into an increase in the effective stressin the granular soil and, in our case, this reflects in aretarded collapse of the granular column. On the otherhand, the loose granular sample undergoes compaction
during the first phase of the collapse, which induces apositive excess pore pressure (see Figure 5). This
positive overpressure decreases the effective stress ofthe granular phase and determines a sudden collapse ofthe column (a liquefaction-like behaviour is observed)
with an extremely rapid moving front.
Similarly to Rondon et al. [3], we measured thepressure at the base of the column (see Figure 1). We
display its time evolution in Figure 6 for both the dense

Fig. 4. Final deposits of the collapsed column for different
aspect ratios: (a) AR = 0.5 and (b) AR = 2. The loose samples
are displayed in the top row (red colour) while the dense
samples are in the bottom row (blue colour).

We observe a good agreement with the experimental
observation. In fact, we obtained a trapezoidal(triangular) deposit for a column of aspect ratio 0.5(AR = 2) when starting from a dense packing fraction.
Instead, we observed a triangular morphology of thedeposit, independently of the column aspect ratio, whenstarting from the loose state.
3.3 Time evolution of the pore pressure
Due to the poromechanical coupling, the evolution of
the pore pressure plays a major role in the collapse
dynamics and one can reasonably assume that it is themain cause of the strong differences observed between

3

E3S Web of Conferences 415, 02025 (2023)
DFHM8

https://doi.org/10.1051/e3sconf/202341502025

The financial support from the Italian Ministry of Education,
University and Research (MIUR) within the framework of theREDREEF Project - PRIN 2017 Call (prot. 2017YPMBWJ)
and the Fondazione Cariverona R&D project are gratefully
acknowledged.

and loose case, referring to a column of aspect ratio
equal to 2, like before.

References
1.

F. Gabrieli, R. Artoni, A. Santomaso, S. Cola,
Collapse of wet granular columns: Experiments
and discrete element simulations, in Particle-Based
Methods II - Fundamentals and Applications, , pp.
562–573 (2011).
2. L. Brezzi, F. Gabrieli, S. Cola, Acta Geotechnica
15(3), 695–714 (2020).
3. L. Rondon, O. Pouliquen, P. Aussillous, Phys.
Fluids 23, 073301 (2011)
4. R. M. Iverson, M. E. Reid, N. R. Iverson, R. G.
LaHusen, M. Logan, J. E. Mann, D. L. Brien,
Science 290, 513-516 (2000)
5. T. Wang, F. Zhang, J. Furtney, B. Damjanac, J.
Rock Mech. Geotech. Eng., 14(3), 1005-1024(2022)
6. B. Chareyre, A. Cortis, E. Catalano, E.
Barthélemy, Transport in Porous Media, 94(2),595-615 (2012)
7. E. Catalano, B. Chareyre, E. Barthélémy, Int J
Numer Anal Methods Geomech, 38(1), 51-71(2014)
8. P. Cundall, O. Strack, Géotechique, 29, 47-65(1979)
9. V. Šmilauer et al. (2021), Yade Documentation
3rd ed. The Yade Project.
doi:10.5281/zenodo.5705394
10. S.C. du Pont, P. Gondret, B. Perrin, M. Rabaud,Phys. Rev. Lett. 90, 044301 (2003)

Fig. 6. Time evolution of the pore pressure at the base of thesample column.

In the first phase of the collapse, we observe a
positive excess pore pressure in the loose case, while anegative pore pressure is registered starting from thedense sample. Moreover, we can observe a longer
transient phase for excess pore pressure dissipation
starting from a dense sample (see Figure 6). We suggestthat this may be related to the larger thickness of thecolumn in most of the stages of this collapse test and to
the lower permeability of the solid skeleton in this case.

4Conclusions

In this work, a numerical study of a column collapse of
granular material in a viscous fluid is presented. ADEM-PFV coupled approach was adopted in order tohighlight the micromechanical aspects which rule thisphenomenon. In analogy with the experimental work of
Rondon et al. [3], the role of the initial solid fraction on
the collapse dynamics was investigated. The numericalresults are qualitatively consistent with experimental
observations highlighting the strong correlation between
the collapse dynamics and the initial packing fraction.
The adopted numerical approach was shown to becapable of capturing the different time evolution of thecollapse as well as the dependency of the morphologyof the deposit on the initial state. Furthermore, the timeevolution of the excess of pore pressure in the basal
section of the column is also well reproduced.
This preliminary study has permitted us to verify thecapability of the DEM - PFV coupled approach toreproduce the collapse of a submerged granular column.
The coupling between the solid skeleton deformationand the pore pressure variation, which is at the base of
the method, well captures the evolution of the porepressure. The DEM-PFV approach seems therefore to bea promising tool for the study of saturated granular
media. However, more studies are needed, in order tounderstand, for instance, the role of the hydraulicconductivity and the contribution of lubrication forces atcontacts.

4