Towards a Micromechanical Understanding of Landslides—Aiming at a Combination of Finite and Discrete Elements with Minimal Number of Degrees of Freedom

In this paper, we propose a combination of discrete elements for the soil and finite elements for the fluid flow field inside the pore space to simulate the triggering of landslides. We give the details for the implementation of third order finite elements (“P 2 with bubble”) together with polygonal discrete elements, which allows the formulation with a minimal number of degrees of freedom to save computer time and memory. We verify the implementation with several standard problems from computational fluid dynamics, as well as the decay of a granular step in a fluid as test case for complex flow.

flow looks much more plausible. Another possible effect which loosens up dense granular matter is "Reynolds dilatancy", where external stresses (in the case of landslides, that may be the additional weight of water due to strong rains) loosen up the soil homogeneously. Any scenario has to obey the laws of physics for particles and grains, and must be verified accordingly. In this paper, we outline a reliable simulation method for particles in fluids in two dimensions. For three-dimensional simulations, the two-dimensional situation should be understood first. A verification by direct micromechanical simulation is preferable to experiments here, as it allows much better control over the initial state and enables us to examine the exact physical processes occurring anywhere at any time within the geometry.

A Combined Simulation Method of Particles and Fluids
As far as the choice between Eulerian (grid-based) and Lagrangian (particle-based) representations is concerned, for the micromechanical understanding of landslides, we chose a particle-based formulation for the soil, and finite elements for the fluid. A list of all symbols with explanations is given in Table 1.

Choice of the Discrete Element Method
For the discrete element method (DEM), often round particles are used, which cannot form heaps (see Figure 1) without unphysical twists (switching off the rotation or introducing an unphysically high coefficient of rolling friction). Also monodisperse particles, especially with negligible elongation, cannot model large granular slopes correctly, as the particles will allow crystal-like ordering which allows slipping along the "crystal axes" with minimal mechanical resistance.
Usually, the most convenient choice to obtain physical behaviour is a distribution of polygons with varying number of corners which are inscribed into ellipses with a statistical variation of the length of the half axes.

Choice of the Fluid Simulation Method
We refrain from simulating the fluid with a particles-based method, because there are problems with the spatial resolution and shot noise would be introduced even in equilibrium flow. For the complicated pore space between polygonal particles, exact discretization is possible via a trigonal grid, so finite elements are the natural Figure 1. The polygonal particles (left) form stable heaps with fluctuations around a "straight" slope. Round particles need support on the sides to form heaps (center) and form a "rounded" slope nevertheless, or they may simply slide along the crystal axes (black lines) and roll away (right) when they have reached the smooth floor. so they offer no advantages in this case). The arbitrary connections of the pore space also require the use of unstructured grids (arbitrary connectivity between elements), though in practice we will generate the grid in the pore space with a constrained Delaunay triangulation, so the connections are limited. Desirable is a formulation with the minimum necessary degrees of freedom, as it is not possible to multiply lattice points near particle boundaries, which would also multiply the computational effort and the memory requirements. Other aspects of the simulation can also be determined as an optimal choice of minimal computational cost, geometrical convenience and the implementation of the actual physical situation.

Interaction between Discrete-Element Particles
The particles in the discrete element method (DEM) interact via a "hard-particle, soft-contact" interaction (see Figure 2): with the characteristic length ( ) proportional to A ∆ , the area change for an integration timestep τ , with the reduced mass of the colliding particles ( ) there is a cutoff so that da el F F < . The solid friction is computed according to Krengel and Matuttis [2], where the static many body friction problem is dealt with "numerically exact" as unilateral constraint problem. The forces between the polygonal DEM-particles are surface forces and are therefore not central, so any torques acting on 1 2 , C C are computed with force arms 1 2 , l l and the sum of the forces (More details about interaction laws can be found in Matuttis and Chen [3]). for the velocities must be one order higher than for the pressures, as the NSE also contains the velocities in one higher order than the pressures [5]. Accordingly, the Taylor-Hood element (P 2 P 1 ) with second order shape functions P 2 for the velocities and first order shape functions P 1 for the pressures had been used in previous simulations [6], but caused issues in case of high particle densities combined with higher flow velocities [7]. Second order elements (polynomials) model a flow field with constant curvature, while for complex flow between granular particles in close proximity (<particle diameter) which move at different or even opposing velocities, the curvature of the flow field can vary. To take this into account, at least third order elements are necessary (see Figure 4). Instead increasing the number of elements between particles has not turned out to be practicable: Using e.g. two elements with different constant curvature side by side does not provide a more plausible flow field or increases the stability of the simulation [7]. Using 4 (or 5) elements would increase the number of mesh points Figure 3. Different triangular finite elements, defined by their polynomial order P. The P 2 + element is a combination of the second order P 2 element and the bubble function (at node J) of the third order P 3 element. Journal of Applied Mathematics and Physics in two dimensions by a factor of 16 (or 25), which means an increase of the computational effort by over one order of magnitude, which is prohibitive. For all of these reasons, in this work, we will be exploring the possibilities of the P 2 + P 1 element, which adds a third order bubble to the second order shape functions ( Figure 3) and thus provides a linear variation in curvature without the need for a full third order P 3 element, while still being LBB stable [4]. While the P 3 looks attractive in theory, for our purposes it has no advantages over the 2 P + element, as we are mostly concerned adapting to the curvature of the flow field.

FEM with Cubic Bubble
On the contrary, having to deal with 10 instead of 7 nodes per element would increase the computational effort and memory consumption by over 40 percent without leading to foreseeable advantages.

Shape Functions
The shape functions are handled in a barycentric coordinate system (local to each element) and must be derived from the Cartesian coordinates of the elements in the grid. Treating every single element separately, we apply a coordinate transfor- , , , x y x y x y of the corners of the triangular element to get the coefficients of its barycentric coordinate functions.

This works easiest via left division (MATLAB's backslash-operator)
with the resulting coordinate functions also acting as the first order (affine) shape functions ψ … so that we have forming the basis of the first order (P 1 ) finite element. Combined with scaling , they act as the element function which represents the pressure p field over our element. Combinations of the affine shape functions ψ … from Equation (4) form the second order (quadratic) shape functions φ … . There are six in total with and their respective permutations, forming the basis of the P 2 finite element. To obtain the 2 P + element, we expand the P 2 basis with the product of all ψ … from Equation (4) 27 the third order (cubic) bubble shape function J θ . J θ in combination with all the φ … functions from Equation (6) and the coefficients ( ) as the element function To treat velocity, a vector quantity, as scalar values, we are splitting it into its Cartesian components, horizontal velocity u and vertical velocity v. Accordingly, Equation (8) represents only the horizontal velocity field and an analogous equ- φ + for the vertical velocity field exists as well. In other words, a shape function can be seen as some sort of unit vector for the field amplitude at its associated node on a single finite element.

Flow Equations
The coefficients ( ) from Equations (5) and (8) With the kinematic viscosity ν η ρ = and the normalised pressure actual p p ρ = obtained by normalizing the dynamic viscosity η and the actual pressure actual p with the density ρ . While rock has a sound velocity which is a multiple of that of water, the sound has to travel through the contacts between sand grains so there is much less material than in bulk rock. Thus, the sound velocity of the granular phase is considerably below that of the fluid and the fluid can therefore be treated as incompressible. The gravitational acceleration g a is also an essential term in the vertical (y direction) Navier-Stokes equation, as the fluid simulation will later run with suspended particles, it is g a which actually leads to the buoyancy force of the fluid acting on the particles. The pressure is calculated from the continuity equation Equations (9)-(11) can now be adapted into their weak form by introducing the arbitrary trial function T for the horizontal velocity, for the vertical velocity and for the pressure. Instead of enforcing a spurious "=0" in the Equations (12)- (14), it is much more meaningful in a numerical implementation (i.e. with rounding and discretization errors) to demand that the right-hand side (the "residual" res) As a rule, a shape function only exists on a single element and is associated to only one node of that element. As a node is usually connected to multiple elements, the combination of all associated shape functions (one for each connected element) is the basis function (in the one-dimensional case this would e.g. be the infamous hat-function) [9]. In our case, the basis function of a node on the edge of a triangle (D, E, F see Figure 3) always consists of two shape functions as the node always connects two elements (except on the boundary). Finally, to approximate the flow field and minimize the residual, an unknown prefactor precedes each test function. When discretizing Equations (12)- (14) for all nodes i, these coefficients must be determined, as they hold the information about the flow state at each individual node. In fact, these are the coefficients from Equations (5) and (8), where (…) stands for the flow variables (u), (v) and (p) at every node i. This results in Equations (12) and (13) and (14)  As the computational effort for the exact spatial quadrature of Equations (12 -14) is prohibitive, we resort to approximation via Gauss quadrature. From the many possible approaches, we opted for the six-point quadrature formula by Gockenbach [10] with a degree of precision ( 4 p deg = ), as any methods with a higher degree of precision produced the same results but with considerably increased CPU-time cost [11]. Using Equation (12) for the horizontal velocity on that node. The notation ( ) is added in the sum of Equation (15). The length of the current timestep is denoted by τ , the index n X − indicates values from previous timesteps. There are certain considerations for the use of the BDF2-solver: While it is A-stable (stiffly stable), so that it can ignore spurious fluctuations from timestep to time-Journal of Applied Mathematics and Physics step (which is an advantage with permanently changing grid geometries), it is also L-stable, which means that it will not damp out instabilities in the flow. This means that BDF2 is not expected to (and, in our experience, does indeed not) delay or suppress the development of vortices in e.g. the flow behind fixed bluff bodies when the flow speed is increased from Stokes flow to Reynolds numbers beyond 10.
The residual of the vertical velocity is calculated analogously to Equations (15) and (16), however the gravity term must be added inside the summation. The pressure residual obtained from Equation (14) is accordingly

Iteration within a Single Timestep
Within a single timestep, Newton-Raphson iteration is used to minimize the residuals. Deriving the residual Equations (15)- (18) for each flow variable gives us the entries for the Jacobian of the Newton-Raphson method [13]. As each equation can be derived by every flow variable (u, v, p) except for Equation (18) where the pressure derivative is trivial, this results in a total of eight different equations for Jacobian entries. Within a single element (with surface e A ), these equations link the flow states u,v and p of two nodes i and j together. They again can be evaluated via the same Gockenbach quadrature method as described above (We omit these equations in this section for brevity and list them in Appendix A of this article). All these entries are combined into the (sparse) Jacobian J via MATLAB's sparse command. As we have largely different scales in the velocities and the size of the finite elements, we use a rescaling of the largest line and column elements of the Jacobian as preconditioning [6].
Using the Jacobian and the residual vector res containing the residuals of all flow states at all nodes, the Newton-Raphson step vector g ∆ for all of these flow states can be calculated as By applying the Newton-Raphson step to the solution vector g 1 g ι ι we complete a full Newton-Raphson iteration ι . Depending on the size of the residual, we now continue with another iteration or advance to the next timestep. As our simulations feature moving granular particles, the geometry will also vary over time. The mesh has to change slightly with every timestep to accommodate the position change of the particles and requires interpolation of the old flow values onto the nodes of the new mesh. Details regarding interpolation of a finite element containing a bubble function are described by Mueller et al. [11].
An overview over the grid generation by regular remeshing via constraint Delaunay triangulation at larger intervals and grid relaxation from timestep to timestep can be found in Ng et al. [14].

Coupling between Fluids and Particles
The simplifying step from three-dimensional to two-dimensional particles means that the DEM-particles are effectively rods in three dimensions. When we assume that these rods would have a rough surface as shown in Figure 5, we end up with a connected fluid space, which can be mimicked in two dimensions by using a larger "shadow" for the particle-particle interaction and a smaller particle "core" which is seen as boundary by the fluid. The width of the core can be adapted to simulate different porosities.
To compute the particle-flow interactions based on hydrodynamic principles, we have to apply no-slip boundary conditions at the particle surfaces. Other approaches, such as immersed boundary conditions [15] violate the no-slip condition and will cause unphysical drag forces between particle and flow. In our approach, the actual force of the fluid on the particle is computed as where the values for the normalized pressure p and velocity vector U are evaluated on the particle boundary Γ from the values computed with the FEM.
Further εζ δ is the Kronecker delta function dependent on the Cartesian directions ε and ζ , n is the normal to Γ and T is the transpose operator. Figure 5. From left to right: Physical three-dimensional particles, three-dimensional rough rods, projection and "averaging" along the z-axis and discrimination into a core region (bold line, dark shading) for the particle-fluid interaction and shadow region (light shading) for the particle-particle interaction so the fluid space in the simulation is connected. Journal of Applied Mathematics and Physics

Verification
To verify the physical accuracy of our simulation code, especially in regard to fluid-solid interaction, we selected three test scenarios. These allow us to confirm that our code is not only providing the accuracy which is expected in the field of computational fluid dynamics, but also that it is capable of handling the specific demands posed by the simulation of granular media and that the continuous change of the mesh does not affect the stability or accuracy.

Wall Correction Factor
We recreated the simulation geometry from Richou et al. [16] to compare our results of the wall correction factor for a fixed two-dimensional circle (respectively a circular cylinder) between two parallel walls. As our simulation is intended to handle polygonal particles, it is impracticable to aim for an implementation of perfectly circular shapes here. Neither curved element outlines nor extremely small discretization would serve any purpose when our main goal is to verify whether the simulation reproduces reasonably well parameters with finite input accuracy. Thus we decided on a regular dodecagonal particle shape as described in the simulations by Ng [6] as a rough approximation for the cylinder.
In the simulations by Richou et al., the main parameters are the Reynolds num- with fluid density F ρ , dynamic viscosity F η and maximum inflow velocity 0 U and the aspect ratio between cylinder size and wall distance b, which both depend heavily on the cylinder radius r. As can be seen in Figure 6, the radius of a regular dodecagon has visible variation which will lead to different results for Equations (22)  (see Figure 6). As buoyancy was not computed in Richou et al. and the particle will not be moving (only the drag forces acting on it will be calculated), gravity was turned off for these simulations.
We have simulated all cases at two different mesh resolutions, first with a thin layer ( max 0.1 r ⋅ ) of mesh elements around the particle and second with mesh elements of roughly the particle edge length (see Figure 7). As Table 2 Table 1.) around the particle. Journal of Applied Mathematics and Physics often is ≥100 elements [15]. This means that for our granular simulations in fluids, the simulation results will not be distorted by a large grid size. The particle shape (which in general will follow a polydisperse size and shape distribution) will be the dominant effect. These results prove the physicality of the microscopic model we use for the fluid-solid interaction as the data by Richou et al.

Cavity Flow
For further verification, we recreated the simulation geometry from Kawahara [18]  Our results for flow field and pressure distribution shown in Figure 8 match the results given in [18] for P 2 b-P 1 (different notation to 2 1 P P + ) to the point, where they are virtually indistinguishable. This is an indicator that our simulation code is not only qualitatively correct, but also quantitatively reliable and able produce results which are expected in the field of computational fluid dynamics.

Granular Step
To prove the capabilities of our simulation code concerning the challenges posed by granular media immersed in fluids, we created a simulation of a collapsing granular step inspired by the experiments performed by Rondon et al. [19]. The findings by Rondon  . To form the granular step, the rightmost column of particles is fixed in place while the remainder of the step is filled by particles falling in via a dry DEM simulation.
When the particles have settled, the FEM part of the simulation is added and the movement of the rightmost particles is unlocked.
As seen in Figure 9, the simulation is able to compute the collapse without Figure 8. Flow field (left) and pressure (right) inside the cavity. These results match the ones generated by Kawahara [18]. issue, while the mesh size in the vicinity of the particles is for the most part roughly the same as the particle edge size (≈1.5 mm). To avoid additional computational effort and unnecessary degrees of freedom, the mesh is coarsened in the region without particles as seen in Figure 10. This grading of the mesh will not have any adverse effects: On one hand, it will retain the original flow field (which would be modified by "cutting off" the fluid domain close to the particles) while on the other hand, the discretization error introduced by the larger elements is compensated by the vanishing of the flow velocity towards the boundaries (Reference runs for a single sinking particle have even shown that the use of a graded mesh gives more accurate results with respect to the symmetry of the trajectory, compared to a fine mesh, because the rounding error due to the conditioning of the Jacobian is reduced [20]).
The particles in our current simulation are still about one order of magnitude bigger than the ones used by Rondon et al., which means we are dealing with a simpler geometry. However this is ideal as a test scenario without sacrificing realism. Now that our simulation code has proven its capabilities the next step is to gradually reduce the particle size (thereby increasing the particle count) to create create a physically meaningful two-dimensional equivalent of the experiments by Rondon et al. and investigate the effect of the packing density on the decay process and the geometry of the end configuration (similar to the simulations by Ng [21]).

Summary and Conclusions
We have presented a simulation of DEM-particles in Fluids where the particles are coupled to the fluid with the full physical interaction (no slip) from computational fluid dynamics. The accuracy of the time evolution was verified via the wall correction factor. For the parameter range under investigation, no instabilities develop, even when mesh elements bordering the particles are about 1/7 th of the particle size. This is a result of the added third order bubble, which stabilizes the simulation for coarser meshes while at the same time guaranteeing a smoother velocity field even with large elements. As a result, less degrees of freedom Figure 10. Full simulation domain for the granular step. Top and right wall are placed far away from the particles to avoid influences from their boundary conditions. The mesh coarsens towards these empty boundaries to reduce computation times. The not ideal shape of triangles in the right upper and lower corner is an effect of the mesh relaxation and the number of prescribed grid points near the right wall. Nevertheless, there are no negative numerical ramifications for the solution from the triangle shape, as the velocity values at the boundary are set and the neighboring computed values are few and of small magnitude. Journal of Applied Mathematics and Physics are needed than those necessary for lower order elements, and the stability is improved compared do P 2 P 1 -elements [7]. In that respect, our simulation is "optimal", as is will not be possible to reduce the degrees of freedom for the FEM further without jeopardizing the stability. As the computational effort is proportional to the number of degrees of freedom, the speed (for a given number of granular particles in a given volume region) will also be optimal.
Neither the accuracy nor the stability is negatively affected by the movement of the particles or the resulting subtle distortions of the mesh in every timestep; the particles come to a physical rest with a physical angle of repose and without any residual noise amplitude. With the current code, we can now move towards larger system sizes.
There remain some "do's and don'ts" concerning the combination of remeshing, graded mesh and finite elements with bubble functions but a discussion of this would lead too far here and will be addressed in a further publication.