CFD-DEM Modeling of Dense Sub-Aerial and Submerged Granular Collapses

Sub-aerial (dry) and submerged dense granular collapses are studied by means of a threephase unresolved computational fluid dynamics-discrete element method (CFD-DEM) numerical model. Physical experiments are also performed to provide data for validation and further analysis. Validations show good compatibility between the numerical and experimental results. Collapse mechanism as well as post-collapse morphological parameters, such as granular surface profile and runout distance, are analyzed. The spatiotemporal variation of solid volume fraction is also investigated. The effect granular column aspect ratio is studied and found to be a key factor in granular morphology for both submerged and dry conditions. The volume fraction analysis evolution shows an expansion and re-compaction trend, correlated with the granular movement.


Introduction
Granular collapses are characterized by the rapid movement of granular material, driven by gravity. They are widely encountered in various geophysical and engineering processes, for example, in the submarine and sub-aerial landslides, debris flows, and river back failures. Understanding the behavior of the granular collapses is crucial for the prediction and mitigation of rated natural disasters with destructive consequences. Nevertheless, granular collapses present very complex mechanics as they often include all regimes of granular flows, i.e., quasi-static, dense flow, and kinetic (or gaseous) regimes in the same setting [1]. The situation is still more complex when the granular collapses are submerged. The interaction with an ambient viscous and dense fluid like water, in presence of the free surface, creates a complex 3-phase system. The fluid inertia, viscosity, and pore pressure can dramatically affect the granular behaviour [1][2][3][4].
There are a large number of experimental [1][2][3][5][6][7][8][9][10] and numerical [1,4,9,[11][12][13][14][15]] studies on gravity-driven granular flows, particularly on granular collapses. While most of the experimental researches focused on the dry (sub-aerial) collapses, a few studies on the submerged condition have been reported (e.g., [1,9]). Earlier experimental studies of Lajeunesse et al. [5], Balmforth and Kerswell [6], Lube et al. [8] investigated the factors influencing the collapse mechanism, runout distance, and the final shape of the granular assembly in the dry condition. They showed that the morphology of the granular mass is largely dependent on the initial aspect ratio of the granular column. Experiments of Rondon et al. [9] showed that the impact of initial volume fraction is more important than the column aspect ratio in case of submerged granular collapses. Granular movement and shear can change volume (and volume faction), which can impact the granular behaviors, especially in the presence of pore fluid in submerged cases [11]. However, this variation volume fraction is not well understood as it has been challenging to quantify in the past researches.
Numerical models have been effective tools for the study of granular collapses in the past. The numerical methods used have been either based on the continuum description [1,[11][12][13]16] or discrete description [4,14,15,17,18] of granular material. The continuum-based models treat the granular material as a body of a complex fluid whose behavior is predicted using a rheological model such as Herschel-Bulkley and µ(I) rheology. Although these models are computationally affordable for large-scale problems, the uncertainties associated with the rheological models may cause inaccuracies, especially in submerged situations. Discrete description, whose most famous example is the discrete element method (DEM) [19], models the individual grains based on Newton's laws. These methods have been widely used and proven their accuracy for dry collapses (e.g., in [17,18]). DEM can be applied to the submerged granular flow by coupling DEM with Lagrangian fluid models or Eulerian computation fluid dynamics (CFD) techniques. However, the application to the submerged granular collapses is limited to only a few studies (e.g., [4,14,15]).
CFD and DEM coupling techniques for granular flows have been either resolved (finegrid) [20,21] or unresolved (coarse-grid) [22][23][24][25]. In the resolved CFD-DEM coupling, the CFD cells are much smaller than the granular particles allowing detailed flow simulation around each particle. Nevertheless, it is computationally demanding. In unresolved CFD-DEM coupling, the CFD computational cell is larger than the granular particles, and the fluid forces are locally averaged over particles. This makes the unresolved methods more computationally affordable. Considering that the granular collapses problems often involve free-surface flows, esp. when collapse-induced wave is important, the CFD-DEM algorithms must deal with the three-phase (gas-liquid-solid) system. A few recent studies (e.g., [24,25]) developed unresolved CFD-DEM methods associated with a 3-phase volume of fluid (VOF) technique to deal with the submerged granular flows with the free surface. Their capability to predict granular collapse remain unassessed.
The present work studies dense sub-aerial (dry) and submerged granular collapses with different aspect ratios using an unresolved CFD-DEM model with three-phase VOF. The numerical modeling is associated with complementary experimentation for validation of the numerical model and further analysis of the mechanisms and processes. In the developed model, the CFD-DEM coupling is based on the OpenFOAM open-source CFD libraries and LIGGGHTS [26] open-source DEM software, coupled through the framework of CFDEMproject [27]. In order to capture the free surface fluid, the volume of fluid (VOF) is integrated into the CFDEM coupling using one of OpenFoam's solvers for modeling multiphase flows. The three-phase VOF is based on the solver cfdemSolverVOF of Jing et al. [24]. The impact of aspect ratio and presence of ambient fluid on the collapse mechanics, morphology, and running distance are investigated. Special attention will be given to the spatial-temporal variation of volume fraction after the collapse, which has not yet been well understood.

Numerical Method
We consider the submerged granular flows with the free surface as a three-phase system that consists of one discrete phase (granular particles) and two continua (fluid) phases of a liquid (water) and a gas (air). The numerical method of this study is based on the coupling of a CFD method, solving the volume-averaged Navier-Stokes equation for the continuum phase (i.e., the pore and ambient fluids), and a DEM method, solving Newton's laws of motion for the discrete phase. Here the coupling of these two methods is based on an unresolved CFD-DEM coupling to reduce the computational time. In this section, these two numerical methods and their coupling is described.

Fluid Governing Equations
The equations governing the unsteady incompressible immiscible two-fluid (waterair) system are the volume-averaged Navier-Stokes equations. A volume of fluid (VOF) method is used to capture the interface of two fluids. In a three-phase system and assuming the unresolved CFD-DEM approach, each CFD computational cell can be occupied by liquid, gas and solid particles as shown in Figure 1. Disregarding the position of particles in each cell, a void fraction for each fluid cell is calculated using a void fraction scheme [24].
Then, α f the void fraction or porosity can be written as: where V pi is the volume of solid particles in each cell, V c is the total volume of each computational cell, and V 1 and V 2 are the volumes of two fluids (water and air, respectively).
We can define V f = V 1 + V 2 as the void space (filled with fluids). In a three-phase system, α 1 = V 1 /V f is the volume fraction of the water over the volume of void. Therefore, α f α 1 = V 1 /V c indicates the local average phase fraction using the transport equation in the VOF method. Therefore, the momentum, and continuity are given as: where u f is the average velocity of the fluid in each cell.
is the stress tensor, f s is surface tension force, F p f is the fluid-particle interaction force. µ f = α 1 µ 1 + (1 − α 1 )µ 2 and ρ f = α 1 ρ 1 + (1 − α 1 )ρ 2 are the volume-weighted viscosity and density of two fluids, respectively. The VOF transport equation, with an artificial compression term, which can lead to sharper interface resolution [28], is given by: where u r = u 1 − u 2 is the relative velocity (also called compression velocity) of two fluids. The fluid governing equations are solved in using finite volume method, with a PISO (pressure implicit with splitting of operators) pressure-velocity coupling algorithm. The VOF equation is solved with the flux-corrected transport (FCT) method [29].

Granular Phase Governing Equations
The motion of discrete phase (granular particles) is described by the DEM method. It solves Newton's second law motion, where the translational and rotational motion of a particle i with mass m i and the moment of inertia I i are governed by [30,31]: where u i and ω i are the translational and angular velocity of particle i respectively, F c ij is the contact force between particles i and j, F f i is the fluid force on the particle i. M ij is the torque from particle j to particle i due to collision.
The contact force F c ij is the core of the DEM method and is consists of two components in the normal (F n ij ) and tangential (F t ij ) directions. It is calculated using the widely used softsphere approach, which assumes a small overlap between contacting rigid particles. Here we use a Hertzian model, which accounts for non-linear elastic contact and damping as: where, δ is particles' overlap, K is elastic stiffness coefficients, γ damping coefficient, and V ij relative velocity. Indices n and t denote the normal and tangential components, respectively. The stiffness and damping coefficients are calculated based on the material properties [31]. The tangential force is truncated to satisfy Coulomb's friction law (F t ≤ µ g F n , µ g being friction coefficient). The interaction forces between the fluid and particles include drag, lift, buoyancy, and the forces due to acceleration of the particle with respect to the fluid (i.e., virtual mass effect and the Basset force) [30]. However, in many multiphase systems, the main and dominant forces are the drag and buoyancy forces, and the other forces are usually considered insignificant [32]. In this study, the drag force and the buoyancy force are considered as the main interaction forces. The buoyancy force is given by: where d p is granular particle size. Here we use the Di Felice drag model [33], which has been used in several of past studies (e.g., [24,34]. In this model the drag force F d is given by: and C d and χ are the drag coefficient and corrective factor, respectively, and are defined as a function the particle Reynolds number. Figure 2 summarized the CFD and DEM algorithms and their coupling. To save computational time, different time resolutions are considered for the CFD and DEM calculations. Given that the DEM usually requires a smaller time step than CFD, the data exchange between the CFD and DEM is only done every n (= ∆t CFD /∆t DEM ) DEM time steps. The computation starts with the DEM solver, detecting the particle positions, contacts, and overlaps. The DEM solves the equation of motion and updates the position and velocity of particles based on the contact, fluid, and gravity forces. Every n DEM step, the particle data are passed to the CFD solver. The updated position of the particles is used to calculate the local volume fraction (or porosity) in each CFD cell, and the momentum exchange information is passed to the CFD side. In the CFD solver, the VOF method is used to capture the free surface and update the volume fractions. Then the pressure and velocity fields are calculated using the PISO algorithm. Fluid forces acting on particles are then calculated and passed to the DEM solver. The model implementation is based on the OpenFOAM and LIGGGHTS, coupled through the CFDEMproject [27], while the three-phase VOF is based on the cfdemSolverVOF solver [24].

Experimentation
Experiments are performed on the dry (sub-aerial) and submerged granular collapses to validate and evaluate the developed numerical model. As shown in Figure 3, the setup configuration includes a horizontal-bed rectangular tank with the length, height and width of l t = 0.30 m, h t = 0.15 m, and w t = 0.1 m respectively. Using a removable gate, the tank is partially filled with granular materials to form a rectangular heap of length l 0 , and height h 0 . In the submerged cases the tank is filled with an ambient fluid (water) of height h f to completely submerge the granular materials. The gate is removed with a speed of around v g = 1 m/s, which allows the near-instantaneous release of granular material. A high-speed camera is used to capture the evolution of granular assembly after release.
The initial aspect ratio of the granular column, a = h 0 /h 0 , is known to have a key role in the collapse mechanism, spreading distance, and the final shape of the granular assembly [3,5,6]. Therefore, we consider two aspect ratios, one below and one above the critical value of 0.7 (suggested by [5]), resulting in a narrow collapse and wide collapse. This results in total of four test cases (D1, D2, S1, and S2) as summarized in Table 1. The test cases have a similar total volume (with around 1% volume difference).   For all cases the granular material is glass-beads with diameter of d = 0.0008 m, grain density of ρ g = 2500 kg/m 3 , a friction coefficient µ g = 0.4 (measured friction angle of 22 • ), and initial solid volume fraction φ 0 = 0.61 ± 0.005, i.e., a dense packing created through tap induced compaction. The glass beads have been colored in a way to form some initially vertical strips to better show the internal deformations. The material has Young's modulus of 60 GPa, restitution coefficient of 0.9 and a Poisson ratio of 0.45. The ambient fluid (for the submerged cases) is water with density of ρ f = 1000 kg/m 3 and dynamic viscosity of µ f = 0.001 Pa.s. Considering the geometrical and material properties, the dry cases (D1 and D2) cases are within the free-fall regime, and the submerged cases (S1 and S2) are within the grain-inertia regime of granular flow.

Dry Granular Collapse
To evaluate the capability of the DEM model, in the absence of the CFD model, it is applied to the dry granular collapse cases D1 and D2. The initial granular configuration is built by feeding the spherical DEM particles from the top into the space bounded by the tank and the gate until the desired packing is achieved. To create the vertical strips (similar to the experiments) the grains of different colors are fed to different regions created by several virtual vertical separators. The material properties are assigned as of the experiments. The gate is lifted with the same speed as experiments allowing the collapse and release of granular material. The results are normalized using the characteristic length of l 0 and h 0 and characteristic time of h 0 /g (the free-fall time scale). Figures 4 and 5 show the snapshots of the experimental and numerical granular assembly as well as velocity field in different dimensionless times ([T] = t/ h 0 /g). They show very similar experimental and numerical configurations. One can differentiate the yielded and non-yielded (immobile) regions from the granular configuration and the velocity field. For case D1 (aspect ratio of a = 0.5), the failure starts with the collapse of the front face. An almost trapezoidal non-yielded region (with a concave front and a flat top) forms at the lower-left corner and grows with time as the yielded region vanishes. Case D2 (with the larger aspect ratio of a = 1.46) has a two-stage process. In the beginning, due to the large height, vertical acceleration is dominant so the failure starts by collapsing the upper part of the granular assembly. Compare to case D1, here the not-yielded region (at the lower-left corner) is smaller and has a triangular shape. Once the overall height of granular assembly is decreased, the horizontal acceleration becomes dominant so the front avalanches and are propagated along the channel. These findings agree with the literature (e.g., [1,5]     To quantify the validations, the numerical and experimental (digitized from imagery data) surface profiles of the granular assembly at different times are compared for two aspect ratios in Figure 6. The figure shows a good agreement between simulations and experiments. The granular evolution is also quantified and validated using the dimensionless runout distances (for different aspect ratios) in Figure 7). We have normalized the results using l 0 (as it has been common in past literature, e.g., in [5,6]), and using h 0 as it might be a more important factor impact the runout. One can distinguish acceleration, steady flow, and deceleration regions in the runout curve as it has also been documented in the literature [1,5,6]. The lower aspect ratio results in a much shorter, and faster-evolving running distance. The numerical and experimental results are in good agreement showing the accuracy of the DEM model of this study. With the normalization of the results using h 0 , the runout distance for the the lower and higher aspect ratios, approach each other. That shows the initial column height might be better normalization factor here. Note that, in the cases we investigate, the total volume is content so change in aspect ratio translates to the change in the height as well.

Dropping Particle in Water
Now that the DEM model is tested and validated, we test and verify the CFD-DEM coupling for a simple widely-used three-phase (air-water-granular) case before applying it to the more complex submerged granular collapse. In this test case, One particle with a radius of 0.1 mm and an initial velocity of zero falls through a free surface (from the height of 1 mm above) and settles inside the fluid (water). The computational domain has a length of 4 mm, a width of 4 mm, and a height of 8 mm. The particle density and air density are 2500 kg/m 3 and 1 kg/m 3 , respectively. The box has partially filled the water with a height of 4 mm and a density of 1000 kg/m 3 . The dynamic viscosity of air and water is considered 10 −5 Pa.s and 10 −3 Pa.s, respectively.
A uniform CFD mesh with a cell size equal to 2.5 times the particle's size is chosen. The simulation results are compared with the analytical solution, based on Stokes's law, taken from [35]. Figure 8 shows a comparison of the particle velocity and the vertical position predicted by the developed CFD-DEM solver and the analytical solution. Although the particle motion in the simulation is slightly delayed compared to the analytical solution (which might be due to the drag model), in general, the simulation results have a reasonable agreement with those from the analytical solution.

Submerged Granular Collapse
Now that the DEM model and its coupling with the CFD model is tested, we apply and evaluate the model of this study for the collapse of submerged granular columns with aspect ratios of a = 0.5 and a = 1.46 (i.e., S1 and S2 cases in Table 1). A uniform mesh with a cell size of 0.0025 m (giving a mesh size to particle size ratio of 3.125) is chosen. The gate is defined in the DEM side only in the simulation. The DEM parameters are the same as the dry cases (D1 and D2). Figures 9 and 10 show the time evolution of numerical and experimental granular configurations, as well as simulated velocity fields two aspect ratios of a = 0.5 and a = 1.46, respectively. As the CFD side of the model does not consider the gate motion, the initial particle suspension and water surfaces disturbance in the experiments can not be observed in the numerical results. A granular failure and movement mechanisms, similar to that of dry cases are observed. Nevertheless, due to the presence of the ambient fluid, the process is much slower, the velocity is smaller, and the front is thicker. Also, the strong counter-clockwise vortex, which is initiated from the top right corner of the granular column and propagates downstream, creates a waveform on the granular surface that is not observed in dry cases. Like the dry case, the column aspect ratio seems to have a key role in the morphology of granular assembly. While for the smaller aspect ratio of a = 0.5 a trapezoidal final deposit (with a flat top) is observed, for the larger aspect ratio of a = 1.46, the final deposit has a triangular shape with a concave top (similar to the dry case). Upon the collapse of granular material the nearby ambient fluid is dragged down and a wave is generated on water surface and propagates downstream and returns after hitting the downstream wall. The higher column aspect ratio cause a higher amplitude wave. In addition to the initial collapse, the wave seems to be affected by the granular movement and the vortex generated by the granular front.
As the column of particles deforms, it drags the surrounding fluid downward and generates forward propagating (seaward) waves and backward propagating (landward) waves when fully submerged.  To quantify the validations, the numerical and experimental surface profiles are compared in Figure 11, for the two aspect ratios, showing relatively good agreement. The discrepancies can be due to the initial gate movement effect on the fluid phase (causing flow disturbance and consequently particle suspension), which is not seen in the numerical model. The evolution of the granular motion is also quantified and validated by plotting the dimensionless runout distances in Figure 12). Similar to the dry cases, three distinctive regimes of acceleration, steady flow, and deceleration are observed, however, the time scale is much larger (close to twice). The numerical and experimental values are in good agreement. With the normalization of the results using h 0 , the runout distance for the tall and wide aspect ratios, get very close.

Volume Fraction Variations
One of the most important factors in the behavior of the material in granular collapses is the particle volume fraction. While the role of initial volume fraction has been well studied in the past [9,11,36], the temporal and spatial variations of volume fraction during the collapse is not well understood, especially in the presence of ambient fluids. The experimental works have challenge estimating the local evolution of volume fraction [7] and the continuum-based numerical models, e.g., [1,[11][12][13], have problems with the accuracy of estimations due to uncertainty with the rheological models. The CFD-DEM model of this study provides a unique opportunity to address this knowledge gap. The average volume fraction of the granular assembly can be estimated by calculating the total volume below the granular surface. The local volume fraction for the submerged cases is directly given by the 3-phase VOF method. Figure 13 provides the evolution of the average volume fraction, normalized using the initial volume fraction, for the dry and submerged cases. Generally, upon release of granular material, a decrease in volume fraction (dilatancy and expansion) is observed. This decrease continues until a minimum value, then is followed by a gradual increase (except for the case S1) or re-compaction. The volume fraction approaches (but never reaches) the initial value. Such behavior correspond to the regime of granular movement (see Figures 7 and 12). Rapid decrease, slow decrease, and gradual increase in volume fraction correspond to the acceleration (and rapid movement), steady flow, and deceleration (and slow movement) of granular material, respectively. For the dry cases, a higher aspect ratio (i.e., case D2) results in a higher variation volume fraction (up to 4%). This is due to stronger acceleration and deceleration, resulting in more expansion and re-compaction of material. For the submerged cases, the behavior is more complex due to the presence of ambient fluid and the material suspension. A similar expansion and re-compaction trend is observed, however, the processes are much longer. The higher aspect ratio (i.e., case S2) results in a higher expansion rate, but eventually, its volume fraction reaches that of the lower aspect ratio.
The local variations of normalized volume fraction for the submerged cases are shown in Figure 14. The values are based on cell volume fraction given by the 3-phase VOF calculations and are only available for the submerged simulations. The dramatic decrease of volume fraction near the granular surface is due to the presence of the water-granular interface, and is not related to granular expansion. As expected, the immobile granular region (trapezoidal in lower aspect ratio and triangular in higher aspect ratio) keeps its initial volume fraction, and the mobile region experiences dilatancy and expansion. The lowest volume fraction (of around 94% of initial volume fraction) is observed in early stages (e.g., T = 2) along the yield curve, where the shear is maximum. Later, the material movement slows down, the volume fraction in the yielded region gradually increases and approaches the initial value. Both cases (S1 and S2) show a similar variation, except that the re-compaction takes longer in case of S2 (with the higher aspect ratio).  Figure 14. Normalized volume fraction field for the submerged cases S1 (top) and S2 (bottom).

Conclusions
An unresolved CFD-DEM model, with a 3-phase VOF technique, based on the opensources platforms Open-Foam, LIGGGHTS, and CFDEMproject, was developed and implemented to study the dense granular collapses in both dry and submerged conditions. The numerical modeling was associate with experimentation, providing data for validation and further analysis. Collapse mechanism and post-collapse morphological parameters, such as granular surface profile and runout distance, were analyzed. Validations showed good compatibility of numerical and experimental results.
The effect of the initial granular column aspect ratio on collapse mechanism and the morphological evolution was significant in both dry and submerged cases. The higher aspect produced a faster running granular front. The two aspect ratios produced different final shapes of granular assembly (trapezoidal for a = 0.50 versus triangular for a = 1.46); however, the difference in final front position was insignificant for the submerged condition.
Investigation of spatiotemporal variation of solid volume fraction showed an expansion and re-compaction trend in the yielded (mobile) region of the granular mass. This trend took longer for the submerged cases. The final volume fraction approached the initial value; however, in non of the cases, the granular mass could not regain its initial volume fraction. The higher aspect ratio resulted in a larger evolution of volume fraction (more expansion and re-compaction) in both submerged and dry conditions. The current study had some limitations that can be address in the future studies. This study did not consider the role of initial compaction (volume fraction) and particle size variations. It was also limited two column aspect ratios. In the model we ignored some the forces of fluid whose effects should be investigated in the future studies.