Micromechanical simulations of the collapse of a submerged granular column

. In this work, we present a numerical study of a granular collapse in a viscous fluid adopting a DEM-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.


Introduction
Wet granular media are frequently encountered in geophysical problems (e.g. debris flows, landslides). As known from classical geomechanics the interaction between the solid and fluid phases plays a fundamental role in the mechanical behaviour of the medium which is mainly controlled by the saturation degree [1] and the fine content [2]. For fully saturated soils, it entails a stress exchange between the phases, which is at the base of 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 to abrupt variations of the dynamics of the collapse in the presence of interstitial fluid. Nevertheless, information at the microscale can hardly be obtained from experimental tests. In this perspective, discrete numerical models represent a fundamental tool and several coupled approaches have been proposed in the technical literature (e.g. DEM-CFD, DEM-LBM, DEM-SPH, DEM-PFV). [5] In this work, we used a discrete element method (DEM) coupled with a pore network technique (PFV) in order 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 capability of the DEM-PFV approach to reproduce the main features of the considered physical process: (i) the collapse dynamics, (ii) the time-evolution of the pore pressure. From this perspective, we compared the

DEM-PFV coupled model
We briefly describe the adopted numerical framework in what follows. The interested reader can refer to [6] and [7] for further details.
The solid phase is modelled through a soft-sphere approach [8] and we assume having a viscous flow of a Newtonian fluid within the pore space. In order to account 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, we assume the two phases to be incompressible and we do not consider lubrication forces. The domain is discretized in tetrahedral regions by using a regular triangulation, in which the vertices of the tetrahedra are the soil particles centres; we refer to the pore space within a tetrahedron as a pore. Considering two adjacent pores, their exchanged flux is linearly related to the local pressure gradient as a consequence of the Stokesian flow regime assumption. In addition, the volume balance equation, for a pore, permits us to link the fluxes to the variation of the pore volume. Finally, by expressing the pore-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 system of equations which, in a matrix form, reads: where K is the conductivity matrix, P is the pressure vector (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 the discrete field of pore pressure. From this, drag forces can be calculated (linearly dependent on P), and these are introduced into Newton's second law and applied to the solid particles.

Numerical results
In this work, we mainly focus on the role of the initial solid 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 Stokesian flow conditions by using the criterion proposed in [10]. Gravity g = -9.81m/s 2 is applied along the y-direction. In the analysis presented here, we mainly consider the following aspects: (i) the time evolution of the profile, (ii) the morphology of the final deposit, (iii) the time evolution of the pore pressure at the base of the sample. This is analogous to the experimental work of Rondon et al. [3]. It should be underlined that the objective of this study is to show the capability of the model to simulate the main aspects involved in the collapse and not to reproduce in a quantitative manner the experimental results of [3]. We will refer to dimensionless quantities in the following analysis. As a length 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.

Time evolution of the collapse
We analyse the influence of the initial solid fraction on the collapse dynamics at first. In this perspective, we focus on the time evolution of the lateral profile and the position 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 the collapsed column is reported every 0.5 s. Starting from a loose sample (ϕi = 0.56) we observe a very rapid collapse dynamics, in which a sudden reduction of the initial height is followed by an almost instantaneous spreading of the column. The front travels very fast and it reaches its final position after a brief time. Conversely, when we start from a dense state of the granular column (ϕi = 0.62), we observe a much slower spreading of the column with a progressive mobilization of the granular material. The collapse phenomenon starts at the top of the column with the triggering of local failures, then it progressively evolves mobilizing the lower layers. We note that the front moves much slower than in the loose case (compare top and bottom profiles in Figure 2) and that the final runout distance is significantly shorter in the dense case. The time evolution 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 a loose initial sample.

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

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 the main cause of the strong differences observed between a loose and a dense initial configuration. In the experimental study of Rondon et al. [3], the basal pore pressure was observed to be strongly correlated to the initial packing fraction of the column.
With this coupled particle-based approach we are able to record the evolution of the excess of pore pressure during the entire collapse and evaluate its distribution in the deformable mass. In Figure 5 we can compare its evolution in the dense and loose sample for the reference case 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 this dilation causes a small suction (pore pressure lower than the hydrostatic one) of the fluid phases from the outer region (i.e., p < 0). Accordingly to Terzaghi's stress principle, the development of a negative excessive pore pressure translates into an increase in the effective stress in the granular soil and, in our case, this reflects in a retarded collapse of the granular column. On the other hand, the loose granular sample undergoes compaction during the first phase of the collapse, which induces a positive excess pore pressure (see Figure 5). This positive overpressure decreases the effective stress of the granular phase and determines a sudden collapse of the column (a liquefaction-like behaviour is observed) with an extremely rapid moving front.
Similarly to Rondon et al. [3], we measured the pressure at the base of the column (see Figure 1). We display its time evolution in Figure 6 for both the dense and loose case, referring to a column of aspect ratio equal to 2, like before. In the first phase of the collapse, we observe a positive excess pore pressure in the loose case, while a negative pore pressure is registered starting from the dense sample. Moreover, we can observe a longer transient phase for excess pore pressure dissipation starting from a dense sample (see Figure 6). We suggest that this may be related to the larger thickness of the column in most of the stages of this collapse test and to the lower permeability of the solid skeleton in this case.

Conclusions
In this work, a numerical study of a column collapse of granular material in a viscous fluid is presented. A DEM-PFV coupled approach was adopted in order to highlight the micromechanical aspects which rule this phenomenon. 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 numerical results 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 be capable of capturing the different time evolution of the collapse as well as the dependency of the morphology of the deposit on the initial state. Furthermore, the time evolution 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 the capability of the DEM -PFV coupled approach to reproduce the collapse of a submerged granular column. The coupling between the solid skeleton deformation and the pore pressure variation, which is at the base of the method, well captures the evolution of the pore pressure. The DEM-PFV approach seems therefore to be a promising tool for the study of saturated granular media. However, more studies are needed, in order to understand, for instance, the role of the hydraulic conductivity and the contribution of lubrication forces at contacts.