A comparison between DEM and MPM for the modeling of unsteady ﬂ ow

. In order to provide a comprehensive comparison between two current numerical methods employed in the modeling of rock avalanches, the Discrete Element Method (DEM) [3] and the Material Point Method (MPM) [1] were used to simulate the mass propagation along a 45° plane transitioning to an horizontal plane. When using the DEM, a 3D code using tetrahedral elements was used and the ﬂ ow was channelized by means of frictionless walls. For the MPM simulations, a 2D code was developed and plane strain simulations were run. Comparisons were made in terms of run-out distance and energy dissipated. In ﬂ uence of parameters such as initial sample geometry, basal friction coe ﬃ cient and shape of blocks composing the sample was studied.


Introduction
Numerical studies involving DEM and MPM were carried out in order to test the ability and the relevance of two kinds of numerical methods (discrete and continuous methods) to reproduce the kinematics of granular avalanches (propagation distances and amounts of dissipated energies). In the field of rock mechanics, the continuum approaches are generally considered when large volumes (i.e. large number of blocks) are involved. For this study, the MPM was chosen over the continuum-based methods because of its hybrid Lagrangian and Eulerian descriptions, which gives it the ability to manage naturally large deformations. This makes it an ideal method for the modeling of granular avalanches involving large masses. However, the continuity assumption may also limit the domain of validity of continuum approaches, specifically when dealing with rock flows. Additionally, the computational times required by the MPM simulations are significantly low, but the dissipations means and more generally the rheologies have to be given explicitly by macroscopic constitutive laws with, unavoidably, the continuity of the medium in mind.
On the other hand, the DEM requires an accurate definition of the initial condition for the terrain, block shapes and families of discontinuities, together with prohibitive computation times. As a consequence DEM turns out to be more appropriate for medium sized volume (i.e. few hundred blocks). Moreover, the complex rheologies inherent in granular materials, are naturally well captured by the DEM -provided that a minimum set of features are taken into account. In our case, these features are the block shapes and adequate dissipative contact models. The price to pay is however long computational times. e-mail: fabio.gracia@3sr-grenoble.fr The aim of this paper is to compare the two types of methods on the base of numerical simulations that consist, first, in the release of randomly packed bricks within a parallelepiped box for the DEM, and second, in a continuum mass governed by an elasto-plastic law for the MPM. The setup is depicted in Figure 1.
A total of 27 configurations were setup using both numerical methods in order to study the influence of parameters such as initial sample geometry, basal friction coefficient and shape of blocks composing the sample (the latter affecting DEM only).

Numerical setup
For all the configurations, randomly packed blocks (cubes of 0.01 × 0.01 × 0.01 m 3 or cuboid elements of 0.02 × 0.01×0.005 m 3 ) were used for the DEM simulations while a continuum mass governed by an elasto-plastic law was considered for the MPM simulations. The container was positioned on an inclined plane at a given height, which is one of the parameters to be studied. Between the slope and the horizontal area, a smooth transition was set by means of a curvature with 0.3 m as radius. DEM simulations were performed in a three-dimensional configuration; the flow was thus canalized to make the comparison possible with the two-dimensional MPM simulations. The canal width was set to 0.2 m to make the flowing truly three-dimensional while limiting the number of elements. This width was chosen so that the ratio between the width and the size of a DEM particle was at least 10. Two different types of containers with varying aspect ratios were used as can be seen in Figure 1. The plane inclination angle was fixed at 45°, while the height at which the block was placed (which refers to the height of the corner closer to the horizontal plane) had three different values: 0.4, 0.8 and 1.6 m.

MPM and DEM laws
As mentioned above, simulations based on the DEM were run using a 3D code while simulations based on the MPM were run using a 2D code. In order to compare the results coming from both methods in an appropriate manner, considerations such as using plane strain when running the MPM simulations and adjusting the MPM material bulk density had to be taken into account to deal with the same mass of material. A classical Mohr-Coulomb elasto-plastic law was used to describe the granular material mechanical behavior. A constant dilatancy angle, close to zero, was used in order to prevent volume increase of the sample. MPM elasto-plastic parameters can be found in Table 1. Within the mass the energy dissipation W is given by: where V p , σ andε are respectively the volume, the stress and the strain increment held by each material point p.
Note that, the index p is not used here for the total stress and total strain to avoid any confusion with plastic terms. Note also that, here and in the following, a parameter super-scripted with a star ( ) relates to the bulk material. Otherwise the parameter concerns the interface between the flowing material and the slope. For the DEM, the interactions between the blocks are governed by a specific contact law established at each contact point, allowing dissipation in the normal and tangential direction of the contact [2] [4]. In the perpendicular direction of the contact a linear elastic law with two different normal stiffnesses (k n ), depending on whether it is loading or unloading, is used. The two stiffnesses relate by using a restitution parameter e 2 n with a value between 0 and 1 which allows for normal energy dissipation. In the tangential direction, a Coulomb frictional model was used, while using k t as tangential stiffness and μ as friction parameter defined as μ = tan φ with φ being the friction angle. The contact parameters governing the mechanical behavior of the DEM elements are given in Table 2.  Table 3. DEM and MPM contact parameters between mass (or blocks) and boundaries at the bottom The elastic contact law and Coulomb frictional model that were used for the normal and tangential contact between particles in the DEM was used to model the boundary conditions (slope, transition and horizontal plane) in both MPM and DEM. The parameters can be found in Table 3. This external force will be added directly to the particles in DEM, while in MPM it will be included in the body forces term involved in the conservation of linear momentum equation.
Concerning the energy dissipation between the blocks (DEM) and the boundary conditions (DEM and MPM), the overall cumulated dissipation was split into 3 contributions: (1) the cumulated work W n of the normal forces (contacts and collisions) at the base of the flow, (2) the cumulated work W t of the tangential forces (friction) at the base of the flow, and (3) the work W of all internal action in the flowing mass. For the DEM, the dissipation within the mass is simply obtained as the sum of cumulated works of contact forces between the blocks in the normal and tangential directions:

Selected results
The numerical configurations were defined by varying the height of release h, the shape of the blocks and the shape of the container (A, B, C), and the basal friction coefficient μ. These configurations are summarized in Table 4. We compared first ( Figure 2) the kinematics of the flow for DEM and MPM in a particular case (simulation #C1, μ = 0.3). We then compared (Tables 5 and 6) the position of the center of mass X CM (called propagation distance here) once the mass had stopped, as well as the distribution of the mass around it. Results were grouped using the basal friction, which is one of the driving parameters in our study. The distribution of the mass was quantified using the standard deviation σ X of the positions of all blocks   or material points depending on the method. In Figure 2 (the figure at the top corresponds to the DEM simulations, while the one at the bottom refers to the MPM ones) we see that the flow evolution from the beginning through the very end correlates properly. Some minor discrepancies can be seen in the shape of the final deposit, as the MPM final deposit spreads a little further. The few blocks that escaped from the mass (spotted at the front of the DEM deposit) are completely absent in the MPM deposit because of its continuous nature. However, even if in some cases there are outstanding similarities, when taking a closer look to the position of the center of mass and the standard deviation around it, the shape of the final deposits are not exactly the same. In particular, the curvature near the slope transition differs between MPM and DEM simulations, and the scattering of the blocks in front of the deposit is completely missed by the MPM. Apart from that, the main mismatch seems to appear at the very end, when the final heap is forming.
All configurations were first tested with a relatively low friction coefficient at the base: the friction angle was set to 16.7°, which is considerably lower than the internal friction angle of the granular material of 28.6°. The  Table 5 for both MPM and DEM simulations. Close correlations between the two methods can be seen.
In regard to the center of mass, X CM , (see Table 5), the value reached by both methods is very similar in most of the simulations. This first observation is rather surprising when one realizes that the constitutive model used in the MPM is rather simple.
In terms of energy dissipation, Table 5 clearly shows that the driving means of dissipation is the friction at the base for both DEM and MPM simulations. This is an expected result since the internal friction angle (i.e., an input parameter of the constitutive model for the MPM, and a parameter mainly related to the friction between blocks in the DEM) is higher than the basal friction coefficient. Internal dissipation achieves nearly the same values in both MPM and DEM. It is however slightly greater for higher container (cases #A). The cases correlating the least (although not largely different) happen to be the ones with higher release heights (cases #B2 and #C2) for which the constitutive model has a larger influence in the stretching of flowing mass. As for #C1, (shown in Figure 2), results correlated properly for the other simulations (from #A1 to #C3). Nonetheless, when using bricks instead of blocks in the DEM simulations, the spreading (σ X ) of the final deposit was reduced; This reduced spreading is an expected result since the rotation of single elements is greatly dependent on its shape, in this case being hindered by the use of bricks. The continuum approach however is not naturally able to simulate it in the way the DEM does.
The container shape becomes important, as well as the block shape, when analyzing the amount of energy dissipated within the granular mass. When using the square container, the energy dissipated in the bulk is larger. These differences might be attributed to the constitutive model itself or to the parameters used since calibration at a very early stage was performed. In order to give a stronger role to the constitutive model, a second series of simulations were carried out with a higher friction at the base.
The friction coefficient at the base was then set to μ = 0.5 corresponding to an angle of 26.6°, which is still below but close to the internal angle of friction of 28.6°. The Table 6 summarizes the results obtained. At first glance it can be seen that similar results are found when comparing DEM and MPM in terms of propagation distance and spreading of the deposit, regardless of the aspect ratio of the starting box, release height or shape of blocks.
When compared with the results obtained with the low basal friction, μ = 0.3, the propagation distances and spreading of the deposit were found to be smaller for all configurations. The observations made with low basal friction are still valid with a basal friction angle close to the internal friction angle.
In regard to the internal dissipation due to plastic straining, an increase can be seen. It is interesting to notice that the better matches are obtained when using a shallow container and cubic blocks for which the assumption of continuity is somehow more appropriate. In other configurations (in particular with bricks), the difference of internal dissipation is more pronounced because of the simplicity of the constitutive model. Despite this, the total amount of dissipated energy is rebalanced thanks to an increase of the friction dissipation at the base.

Conclusions
The use of the material point method has been evaluated in the case of a flow in the transient regime that involves finite strains. In the comparisons presented, the DEM results were taken as the ground truth. Although an extensive simplification has been considered in the continuous model, the MPM has proven to be of interest in the modeling of large size rock avalanches or landslides. In terms of energy dissipated, with the assumptions made for the constitutive law of the continuum model (Mohr-Coulomb and values of parameters) some discrepancies started to become evident as the basal friction coefficient was increased (which involves shearing within the granular mass), but the results were still related.
In the field of soil (hydro-)mechanics or computational fluid dynamics, many constitutive models exist. The MPM provides a framework that is not more restrictive than all approaches based of the continuum assumption. Accordingly, some events such as debris flows, mud flows or snow avalanches -deemed more complex to model -can be dealt with by implementing the best suited constitutive models. A more ambitious project is to make use of twoscale modeling by replacing the constitutive law by DEM simulations (representative to the material) for each material point.