A partitioned material point method and discrete element method coupling scheme

Mass-movement hazards involving fast and large soil deformation often include huge rocks or other signiﬁcant obstacles increasing tremendously the risks for humans and infrastructures. Therefore, numerical investigations of such disasters are in high economic demand for prediction as well as for the design of countermeasures. Unfortunately, classical numerical approaches are not suitable for such challenging multiphysics problems. For this reason, in this work we explore the combination of the Material Point Method, able to simulate elasto-plastic continuum materials and the Discrete Element Method to accurately calculate the contact forces, in a coupled formulation. We propose a partitioned MPM-DEM coupling scheme, thus the solvers involved are treated as black-box solvers, whereas the communication of the involved sub-systems is shifted to the shared interface. This approach allows to freely choose the best suited solver for each model and to combine the advantages of both physics in a generalized manner. The examples validate the novel coupling scheme and show its applicability for the simulation of large strain ﬂow events interacting with obstacles.


Introduction
Mass-movement hazards involving fast and large soil deformation have increased significantly during the past decades due to climate change and global warming. Those phenomena like debris flow, avalanches and mudflow can cause extensive damage to landscapes and infrastructures. Even more devastating are those phenomena, which carry huge rocks or other significant obstacles, causing high economic loss and often human casualties. Therefore, further assessment on the evolution and effects of such disasters as well as the definition of predictive tools and of countermeasures are in high demand. One of the critical issues for further understanding of such phenomena, as well as the prerequisite for the construction of countermeasures, are numerical simulations.
For the simulation of those large strain flow events, which include topology changes of the material, standard discretization techniques such as the Finite Element Method (FEM) are likely to suffer from mesh entanglement and distortion requiring expensive remeshing additional material points for the representation of each discrete block are introduced so that the contact can be detected through the nodes of the computational background grid. The resulting acceleration at the grid nodes arising from the interaction is projected to the nodes of the discrete element working as body forces, whereas the contact between the blocks is detected using a shrunken point method. The proposed scheme is enhanced by Jiang et al. [37] using a spheropolygon DEM (SDEM) [38] which is coupled with MPM by exchanging contact forces. As the contact detection and force evaluation is unified under the scheme of SDEM, the dependency on the computational background grid is reduced significantly compared to the previous approach, and the influence of different particle shapes can be better preserved. During the calculation, MPM particles which are in the Verlet distance, which is the cut-off distance of the potential contact between two discrete elements, are treated as small SDEM disk particles with a certain radius for contact detection resulting in DEM contact forces, which are applied in the form of an additional boundary force term to the material points.
In this article, we propose a novel partitioned MPM and DEM coupling scheme, which allows us to use the best-suited solution strategy for each sub-problem, whereas the interaction is shifted to the shared interface leading to an exchange of data. Therefore the user is not restricted to a code that includes both participants, and even more important, no special treatment as transforming material points to DEM particles [37] or introducing additional material points at the outline of DEM particles [36] for the contact detection have to be considered. Instead, it is possible to couple any existing DEM and MPM software by creating a suitable interface, leading to an efficient and robust partitioned coupling scheme. For this work the open-source multiphysics software KRATOS [39][40][41] has been used, which is written in C++ and offers a Python interface.
This paper presents the basic theory, introducing the used notation of MPM and DEM in "Material point method (MPM)" and "Discrete element method (DEM)" sections, respectively. Then, in "Partitioned weak MPM-DEM coupling scheme" section, the partitioned coupling scheme of DEM and MPM is presented and finally validated in "Examples and validation" section. Several examples are conducted, starting with a beam impacted by a DEM sphere, which can be solved analytically, and therefore the numerical results of this example constitutes the validation and verification of the proposed coupling scheme. This example is calculated in two-as well as three dimensional case. Additionally by enlarging the beam, a multiple impact scenario can be created and the results can be compared to the literature. Finally, to show the applicability for mass-movement hazards, the numerical results are compared to the experiment conducted by Liu et al. [36] as well as to reference solutions [36,37] from literature.

Governing equations
In a Lagrangian kinematic description of a continuum body B, which occupies a domain in a three-dimensional space E the conservation of mass and linear momentum needs to be satisfied, assuming isothermal settings. In these equations, ρ represents the mass density, b the volume acceleration and σ is the symmetric Cauchy stress tensor. The first and second material derivatives of the displacement u are the velocity and the acceleration, respectively. The problem (1)-(2) is fully defined with the boundary conditions where u is a prescribed displacement on the Dirichlet boundary D and p is a traction vector on the Neumann boundary N with normal n.
Multiplying the momentum Eq. (2) with a virtual displacement field δu and integrating over the domain ⊂ E leads to the weak form of the balance equation: with the virtual strain δe arising from the virtual displacement field δu. Thus Eq. (5) can be expressed as a variation w.r.t. δu as Considering geometric and material nonlinearities, a linearization of the weak form is necessary, and thus the Newton-Raphson method is used to approximate the solution iteratively.

Discretization in time and space
In MPM, the body B is discretized (indicated by superscript h) into a finite number n p of Lagrangian moving particles, also called material points: Those particles, assigned a constant mass m p , represent a finite volume of the body p , and they carry the history-dependent variables and the material information during the calculation process. Additionally a fixed background grid is introduced to approximate the continuous fields and their respective gradients by locally defined basis functions N I and discrete quantities at the respective nodes I of the corresponding background element. Therefore Eq. (6) is rewritten in its discretized version: where R represents the residual force vector which has to vanish, since the virtual displacements are arbitrary. The linearization of this non-linear equation leads to the discretized solution system where K is the tangential stiffness matrix and u h the incremental displacement for which the system is solved. The entries of the tangential stiffness matrix, related to the nodes I and J of the background grid element, as well as the residual force vector components, related to the respective background grid node I, can be obtained by: where D is the constitutive matrix, B I the deformation matrix related to each node I in the element and ∇ x denotes the spatial gradient at the current configuration. The operator specifies the assembly procedure that takes into account the kinematic compatibility, while the factor α in Eq. (10) is defined by the integration scheme which is defined by α = 1/(0.25 t 2 ) applying an implicit Newmark time integration scheme. Notice, that the traction surface integral in Eq. (11) is kept in the weak form, as the boundary discretization is discussed in "Boundary conditions in MPM".
Due to this dual description of material points and computational background grid, data from the material points (marked with underscore p in the equations) has to be extrapolated to the nodes of the computational background grid (marked with underscore I in the equations) and vice versa, resulting in the three phases of a MPM calculation, also illustrated in Fig. 1: II Lagrangian phase: Solution of the discretized governing equations, resulting in a deformation of the background grid. This step is similar to an Updated Lagrangian FEM step, therefore the reader is referred to [42] or [43][44][45] for a detailed description. III Convective phase: Solutions obtained at the nodes are interpolated back to the material points, resulting in an update of the material point's position, and the background grid is reset to its initial position.
The summation to receive the position and acceleration update is performed for all connectivity nodes n n whereas the velocity is updated via a trapezoidal rule where the upper index n indicates the time-stepping and t the time-step size.
The implicit Newmark time integration scheme applied in this paper is an extension of the initially proposed displacement-based formulation presented in [46]. Further details, as well as an extension to incorporate mixed formulation, can be found in [42,47,48].

Boundary conditions in MPM
Boundary conditions in MPM can be applied directly at the nodes of the computational background grid if the boundary of the material domain matches with the background grid. This means that both Neumann and Dirichlet boundary conditions can be imposed in a FEM fashion, complying naturally with the Kronecker delta property. However, as the material points move independently of the grid, the boundaries do not necessarily coincide with the nodes of the background grid leading to the necessity of imposing in-homogeneous or non-conforming boundary conditions for the general case. Several attempts have been made to address this issue, just mentioning among them the Penalty Method [49] for the weak imposition of Dirichlet conditions, and [50] to impose moving Neumann boundary conditions. Focusing on the imposition of point loads for the proposed coupling scheme, boundary particles are introduced [51] to discretize the continuous Neumann interface N by a discrete number of mass-less particles n bp : where A bp indicates the current area of the boundary particle. Hence, the traction surface integral of Eq. (11) can be rewritten as where F S is the respective resulting point load at each introduced boundary particle. The subscript S, which is introduced at this point, indicates the imposition of the point load within the MPM solver. Including this in the MPM calculation procedure means, that in the Initialization phase a search for each boundary particle is performed to define the connectivity before the respective point load is mapped via the corresponding basis functions to the nodes of the background grid.
Then the system is solved in the Lagrangian phase before the kinematic variables at the boundary particle are updated, following the concept of material points. Therefore the boundary particles move according to the deformation of the body discretized by material points. This is an important feature for the proposed partitioned MPM-DEM coupling scheme, as those boundary particles are tracking the interface of the MPM domain.
Special consideration is necessary for background elements, which carry boundary particles but no material points. Therefore, to fulfill the equilibrium condition, the point load values are mapped solely to nodes of the background grid, which are assigned a mass and therefore are connected to the body. This can be achieved by modifying the basis functions N I (x bp ) of the nodes evaluated at the position of the boundary particle x bp : where is usually considered numeric zero whereas the partition of unity is guaranteed by the weighting procedure.

Discrete element method (DEM)
In contrast to MPM, which belongs to the group of continuum-based particle methods, DEM considers the motion and interaction of discrete particles. Since its first mention and derivation in [29], DEM has become increasingly popular and is now used in both industrial applications and science. DEM is a particle method, which discretizes the considered body into discrete particles. Therefore it is often used for granular materials on one hand and interacting discrete obstacles on the other side. Additionally, large discrete events, such as rockfall [30][31][32] can be efficiently handled with the DEM. While arbitrary polyhedral shapes of the particles could be considered [28], contact algorithms applying a sphere are preferable, as the contact calculation is efficient. For the simulation of arbitrarily shaped bodies, the clustering of DEM particles which has been investigated in [52] will be applied. The creation of such clusters has been extensively discussed in [53,54], which provide a free-to-use online tool [55].

Contact detection
As described in [28] this approach results in a reduced computational time when calculating overlaps compared to the contact between two polyhedra. Using spheres or clusters of spheres, only the contact between sphere-sphere, sphere-line, sphere-vertex, and spheresurface need to be considered, which are simple operations in which only the shortest distance and the respective sphere radius are compared. This has been investigated in [33,34], which additionally describes an efficient way of handling various contact partners at the same time, applying the so-called Double Hierarchy Method. A sphere with the center C i and a corresponding radius R i is in contact with an arbitrary geometric object as soon as the shortest distance d i from the surface of the object to the center of the sphere C i is smaller than the radius, that is d i < R i . See [34], on how to calculate d i for different geometric entities, such as vertices, lines, and surfaces.

Contact forces
As soon as contact is detected, the respective contact forces are calculated. A variety of different contact laws can be applied, while a Hertz-Mindlin spring-dashpot (HM+D) [56] model is used in this work.
As Fig. 2 visualizes, the contact between different geometric objects with the HM+D model needs the definition of various DEM parameters: • k n , k t : Normal and tangential spring stiffness.
The proceeding calculation of the contact force, with the help of the DEM parameters mentioned above, is omitted at this point and can be looked up in [33,34,[56][57][58][59].

Integration of motion
After computing the contact forces, the DEM solution process proceeds to the integration of motion, following Newton's second law of motion. While the mass m relates the translational accelerationü to the forces F, the inertia tensor I is used to calculate the moments (torques) T via the rotational accelerationω.
Following [30,34,60] the forces and moments on each of the particles are described as the sum of multiple contributions, With, • F ext i , T ext i : External loads, torques. • F ij : Contact forces between sphere i and sphere j or MPM boundary j as F i, Contact .
• n ij , t ij : Normal and tangential vectors at the contact points between the two respective geometric objects i, j. • r ij c : Connection vector between sphere i and the contact point to the neighbour j.
Different time integration schemes can be applied, ranging from a first-order forward Euler to arbitrarily high-order schemes with an increasing level of complexity. For the simulations in this work, a second-order Velocity-Verlet [28] (central differences) scheme is used to integrate the translational degrees of freedom. In contrast to the classical Verlet method [61] the Velocity-Verlet provides second-order accuracy for both the displacement as well as the velocity. Furthermore, to provide a robust time integration of the rotations [62] proposes a time integration, working with quaternions [63] which is used in the following examples.

Partitioned weak MPM-DEM coupling scheme
The fundamental idea of a partitioned or staggered coupling scheme is that the involved sub-solvers, specifically for the proposed coupling scheme the MPM-and DEM-solver, are treated as black-box solvers and the communication between them, the MPM (indicated by subscript S) domain S with boundary S and the DEM (indicated by subscript P) domain P with boundary P is shifted to their shared interface where the interface transmission conditions need to be satisfied. Enforcing the consistent deformation of the involved sub-systems leads to the kinematic condition where u S (x), u P (x) are the continuous displacement andu S (x),u P (x) the velocity fields of the MPM and DEM partition at the interface respectively. Additionally, the load balance has to be satisfied, leading to the dynamic condition where p S (x) and p P (x) are the traction fields at the respective interface, which are defined with respect to the corresponding outward normal vectors. Discretizing the involved subsystems for the solution process, the kinematic constraint is formulated in the discrete form as where u P , u S are the displacements andu P ,u S are the velocities at the respective nodes of the discretized interface. The matrix H PS is the direct mapping matrix, arising from the applied mapping technique [64,65], as the discretizations of the sub-systems are in general non-matching.
Rewriting the dynamic condition in its discrete form results in where F S are the resulting external forces at the MPM interface nodes and F P the resulting contact forces at the nodes of the DEM partition. In this study a conservative mapping technique [64,66] is applied and therefore the transposed mapping matrix H T PS is used, to map the forces from DEM to MPM, enforcing the conservation of energy at the shared interface.
In the partitioned strategy the two partitions are solved separately and the interface equilibrium is enforced by boundary conditions which are imposed at each interface. Therefore in the DEM partition a Dirichlet condition, called wall condition in the following is introduced at the shared interface which was originally proposed by Santasusana [33] for partitioned coupling with FEM. Depending on the velocity and the displacement of this wall condition, contact forces with the DEM particles can be calculated. During the calculation process this interface is discretized into line segments whereas the respective surface is triangulated in the three-dimensional case. Therefore, solving the DEM partition possibly results in contact between the DEM particle and a vertex, a line segment or a triangular surface, leading to contact forces F i,Contact at the DEM particle i and discrete forces F P at the corresponding nodes of the wall condition.
Thus, expressing the DEM partition as black-box solver results in where u P ,u P represent the discrete nodal displacements and velocities of the wall condition, which are the input for the DEM solver and finally resulting in contact forces F P at the discretized nodes of the interface condition after solving the DEM problem (details of the DEM calculation process are summarized in "Discrete element method (DEM)" section). Applying a Dirichlet-Neumann partitioning, a Neumann boundary condition is introduced in the MPM model at the shared interface. This interface is discretized by boundary particles which can be assigned a point load and are kinematically updated according to the material points, tracking the interface. Therefore, meeting the discrete dynamic interface condition (Eq. 23) the contact forces F P calculated by the DEM solver are mapped to the MPM partition as external forces F S . For the data transfer between the two interfaces the nearest neighbor interpolation technique [64] is applied in this study resulting in a copy operation in most of the cases, as the boundary particles usually are located at the same position as the nodes of the discretized wall condition in DEM to reduce errors arising from the mapping.
Therefore the equation for the MPM black-box solver can be written as (details of the MPM calculation process are summarized in "Material point method (MPM)") taking the forces F S as input. The solution of the MPM domain then results in new displacements u S and velocitiesu S at the boundary particles which are discretizing the shared interface within the MPM partition. Meeting the discrete kinematic constraints (Eq. 22) this interface update consequently leads to updated displacements u P and velocitiesu P as an input for the DEM solver. Thus, sequentially solving Eqs. (25) and (26) which are determined by the interface transmission conditions (Eqs. 22 and 23), result in the Gauss-Seidel communication pattern [67,68], which is schematically visualized in Fig. 3 for the weak coupling scheme. In detail, bringing all aspects together the following workflow results for a complete partitioned weak MPM-DEM coupling scheme: First of all, the models for the MPM and DEM partition are created independently of each other. At the shared interface, which usually coincides with the outline of the MPM body, a wall is defined in the DEM partition whereas boundary particles are introduced in the MPM counterpart. Avoiding errors due to data mapping, it is recommended, to locate the boundary particles in the MPM model at the same positions as the nodes of the discretized DEM wall.
After the problem setup, the calculation procedure starts, which is also summarized in algorithm 1 and visualized in Fig. 4, assuming that both involved solvers advance in time with the same time-step. The DEM partition is solved first, with given displacements u P and velocitiesu P at the nodes of the discretized wall condition, possibly leading to some contact between the DEM particles and the wall condition. The resulting contact forces at the nodes of the discretized wall condition are then transferred by Eq. (23) to the interface of the MPM partition as external forces F S .
In the MPM calculation process, those point loads F S at the boundary particles are treated as non-conforming Neumann conditions, leading to the calculation process, described in "Boundary conditions in MPM". According to the deformation of the material points, also the displacements u S and the velocitiesu S of the boundary particles are updated, leading to a kinematic update of the shared interface. This update is mapped by Eq. (22) to the DEM partition updating the nodal displacements u P and velocitiesu P of the DEM wall condition accordingly.  (2). Then MPM is solved leading to a kinematic update of the MPM interface (3), which is mapped back to DEM interface (4). Due to weak coupling, the steps are repeated for the next time steps Applying a weak coupling scheme also known as an explicit coupling scheme [22], the DEM solver advances in time and solves the DEM partition with the updated wall. For each time step, those steps are repeated until the end of the simulation t end .
It is important to remark that a DEM particle i calculates one resulting contact force F i,Contact , depending on the shortest distance between the interacting object and the center of the particle. This can lead to difficulties during the coupling simulation especially considering solid or granular material in the MPM domain and their interaction with rather large obstacles compared to the background mesh size.
As illustrated in Fig. 5a the DEM particle i calculates a single contact force F i,Contact with the wall condition in each time step and resulting in forces F P at the corresponding nodes of the discretized wall condition. Assuming a two-dimensional case and a DEM particle interacting with a line segment, this would result in two forces at the nodes of this segment. Mapping those forces to the MPM partition leads to point loads at the boundary particles, and usually, as a matching discretization of the interface is used, it would result in two point loads. Depending on the location of the boundary particles within the background mesh and the element size of the mesh, a particular area of the MPM domain is affected by the interface forces within each time step. Especially for granular material or solids within the MPM domain, together with a relatively larger particle size compared to the element size of the computational background grid in the MPM domain, this could lead Calculation of contact forces in DEM and influence of the MPM partition. a One DEM particle calculates one resulting contact force, affecting only a limited number of background grid nodes in the MPM partition. b By clustering the DEM particles, also suited for arbitrarily shaped particles, the contact forces are distributed over the interface, leading to an accurate interface representation to an insufficient representation of the shared interface, possibly resulting in penetration of the MPM material into the DEM particle and thus negatively affecting the results of the coupled simulation.
To resolve this issue, clustering of DEM particles [52] can be applied. Besides the advantage that arbitrarily shaped obstacles can be modeled efficiently by clustering DEM particles, each DEM particle within the cluster can calculate a contact force, leading to a more accurate representation of the interface, depicted in Fig. 5b. Therefore the clustering method should be applied if the shape of the particle should be accurately considered and a detailed interface representation is necessary in the MPM model.

Examples and validation
A series of tests have been conducted to show the application of the proposed partitioned MPM-DEM coupling scheme and its accuracy. First of all, the results are compared to an example with analytical solution from the literature for two-and three-dimensional problems. The second example which is an extension of the first one, represents a more challenging impact scenario. Also for this example a reference solution found in the literature is used for the assessment of the proposed numerical approach. Finally the numerical solution of a granular flow impacting on wooden obstacles is compared to experimental results conducted by Liu [36] to show the application for Soil-Particle interaction.

Single impact of DEM particle on simply supported beam modelled with MPM
For the validation of the proposed algorithm, an academic particle-structure interaction example is chosen, which was first proposed by Timoshenko [69] in 1951 and reviewed in detail by Meijaard [70]. It consists of a simply supported beam, impacted at its center by Fig. 6 Single spanned beam hit at its center by a sphere a sphere, and was also studied by Santasusana [33] to validate the partitioned DEM-FEM coupling scheme. The beam, in [33], was modeled by FEM solid elements, whereas a DEM particle represented the impacting sphere. According to this setup, in the proposed study, the impacting sphere is also represented by a DEM particle, whereas the elastic beam is calculated, using MPM.
The flexible beam, single supported on both sides, has a total length of 15.35 cm and a square cross-section of A = 1 × 1 cm 2 whereas the impacting sphere has a radius of R = 1.0 cm and an imposed velocity ofu 0 = −0.01 cm s in y-direction. The system as well as the material parameters described in [70] are illustrated in Fig. 6.
Two models are created to validate the proposed coupling scheme and the respective boundary conditions for both the two and three-dimensional case. For the discretization of the two-dimensional MPM model, a quadrilateral background grid is chosen, which has three times the height of the beam and a total length of 15.87 cm meshed by 62 × 13 elements. Therefore the resulting element size of the computational background grid is similar to the mesh size chosen by [33] when discretizing the beam by FEM solid elements. To generate the material points, an additional quadrilateral body mesh with 120 × 8 elements is created in the preprocessing step. This mesh is used to initialize the material points at the respective Gauss point positions. For this example 16 particles within each body mesh element are considered, minimising the error arising from particle integration. The shared interface, defined at the top edge of the beam, is divided into 60 line segments for the DEM wall condition, whereas in the MPM model, 61 boundary particles are initialized coinciding with the nodal positions of the DEM wall.
The wall condition in the three-dimensional example, introduced at the top surface of the beam is discretized by 77 × 10 structured triangles. Again, the boundary particles are placed at the same positions as the nodes of the wall, reducing errors arising from the mapping. For the three-dimensional MPM model, a structured hexahedral mesh is considered similar to the two-dimensional case but with a width of 1.52 cm being discretized by seven elements following [33]. For the body mesh the width direction is divided into 10 elements resulting in a hexahedral mesh. Again, 16 material points are initialized at the respective Gauss point positions of the body mesh elements.
For the impacting sphere, a Hertzian contact law is considered, and the coefficient of restitution is chosen to be 1.0, whereas a zero friction coefficient for the interaction Fig. 7 Single spanned beam with a total length of 15.35 cm hit laterally at its center by a sphere with radius R = 1 cm. Numerical results for 2D and 3D simulation in comparison to reference solution [70] between DEM and MPM is considered. For both examples, a time step of t = 5e −8 is chosen.
During the calculation, the displacement of the sphere and the deformation of the beam center is measured. Additionally, the total contact force resulting from the impact is observed during the calculation time. The results are printed in diagram 7 and show excellent agreement with the solution obtained by [70] and [33].
Due to the chosen geometry, this example produces a single impact between the beam and sphere, resulting in an oscillation of the beam, depicted in blue in the diagram 7 which corresponds to the natural frequency of the structure.

Multiple impacts of a DEM particle on simply supported beam modelled with MPM
The configuration of the second validation example is similar to the previous one except for the length of the beam which is increased to 30.70 cm and the radius of the sphere being R = 2.0 cm in this case resulting in multiple impacts of the DEM particle on the simply supported beam. Also for this case, the reference solution can be found in [69], [70] and [33] investigated this example to validate the proposed DEM-FEM coupling scheme, too.
Following [33] the background grid is meshed by 124 × 13 quadrilateral elements considering a total length of 30.96 cm for the grid. The material points are initialized at the Gauss point positions of a structured body mesh with 24 × 8 quadrilateral elements. Also in this example 16 particles within each body mesh element are generated, reducing the error arising from particle integration within the MPM calculation procedure. Again, the top edge of the beam is the shared interface which is divided into 240 line segments for the DEM wall condition while the boundary particles in MPM are placed at the same positions as the nodes of the discretized DEM wall. Single spanned beam with a total length of 30.70 cm hit laterally at its center by a sphere with radius R = 2 cm. Numerical results compared to reference solution [70] and numerical results of a DEM-FEM coupling scheme conducted by [33] The obtained results for the displacement of the beam center and the sphere as well as the resulting contact forces during the simulation are plotted in Fig. 8 in comparison with the reference solution proposed by [70] and [33].
Due to the changed geometry of the system, compared to the first example, the sphere impacts the beam three times leading to the vibration of the beam. Also for this example, the results show a very good agreement with the analytical solution proposed by [70], however, a small difference of the data during the second and third impact can be observed. This effect was also observed by [33] where the beam was modelled by linear solid elements (FEM), which similarly to the MPM model are not correctly capturing the higher vibration modes and therefore leading to the small deviation on the second and third contact events. Therefore, compared to the numerical solutions proposed by [33] a nearly perfect agreement can be observed. The slight difference in the two simulations are inherently arising from the MPM calculation, as the nodes of the computational background grid are not coinciding with the boundary of the considered beam, leading compared to FEM to a cruder approximation of the displacement field and a weak imposition of the contact forces.

Granular flow impacting DEM obstacles
The third validation case considers the interaction of significant strain flow events with obstacles. For this purpose, the experiment of granular material impacting wooden blocks, which was initially conducted by [36] is simulated and the obtained results are compared with the available data from the literature at specific time steps. Furthermore, this experiment was also investigated numerically by Liu et al. [36] for the validation of a monolithic coupled MPM-DEM method as well as by Jiang et al. [37] validating their proposed MPM-SDEM coupling method. In both numerical simulations from the literature the granular flow is modeled by MPM whereas the discrete wooden blocks are represented as DEM [36] or SDEM [37] particles. In the coupling scheme proposed by [36] additional material points at the outline and the center of the discrete wooden blocks are introduced, so that the contact between the discrete objects and the granular flow is simulated by a momentum exchange at the nodes of the computational background grid. Therefore the accuracy of the method hardly depends on the background grid size which is improved by [37], where both the contact detection and the force calculation is unified under the contact scheme of the SDEM method. In this case, however, it is necessary to transform the material points which are within the Verlet radius of a SDEM particle into small SDEM disk particles with a certain, user-defined radius to calculate the contact forces which are applied to the material point as an external boundary force. Both methods are using a monolithic approach for the coupling of MPM and DEM and the respective numerical solutions are compared to the simulation results of the proposed partitioned MPM-DEM coupling scheme.
The initial configuration of the two-dimensional model representing the experiment is depicted in Fig. 9.
consisting of the granular material confined initially to a region of 10 × 20 × 20 cm 3 and after releasing flows down due to gravity. Additionally, three identical wooden blocks are considered, which are placed on top of each other at a distance of 30 cm of the granular material, whereas the block on the bottom is glued to the desk to simulate a foundation of a building. Due to gravity, the granular flow impacts those wooden blocks, pushing the upper two wooden particles to the right side, resulting in an angular velocity, whereas the block on the bottom stays fixed.
For the calculation of this experiment, the bottom of the box containing the material is modelled as a fixed boundary while a slip condition is assumed for the vertical wall. The granular flow is simulated by MPM using Mohr-Coulomb plane strain material law with the respective material parameters summarized in Table 1 taken from [36].
For the discretization of the MPM model, an unstructured triangular background mesh with an element size of 0.5 cm is chosen, whereas for the initialization of the material points a body mesh of structured triangles with half of the element size containing three particles each is selected. Compared to [36] where the background grid size is determined by the contact detection, a finer discretization is chosen for our MPM model, to accurately reproduce the movement of the granular material between the wooden blocks in the run out of the experiment. As the granular material undergoes large deformation, boundary particles are placed with a distance of 0.005 cm around the initial configuration of the granular material to ensure a suitable discretization of the shared interface throughout the simulation. For this example however, those boundary particles are not placed exactly at the outline of the body but are shifted marginally inside the body, to avoid numerical instabilities during the calculation. Considering gravity driven granular material, this modification is necessary, as the boundary particles receive the resulting contact forces from the DEM partition and apply them to the MPM model as point loads. If those point loads are applied to elements which contain only few single material points-which is a general case at the body outline-this could lead to non-physical behaviour of those elements and therefore to large movement of the material points within those elements. Therefore to ensure that the forces are applied to the main material flow, the boundary particles are defined within the first row of the material points. In this specific case, the marginal shift δ is one third of the body mesh size, much smaller than the size of the background elements, and therefore almost negligible for the final solution. In addition to the marginal shift δ, future research could investigate alternative solutions to overcome the numerical instabilities at the interface such as [71,72].
The wooden blocks are simulated by DEM, whereas each block of size 2 × 1.8 × 19.8 cm 3 consists of 8 × 8 spherical particles which are compacted to a cluster to model the squared shape of the blocks. As the block on the bottom is glued to the desk, it can be modeled as a fixed boundary in the simulation, and therefore only block one and two are represented as a DEM cluster particle. The discretization of the boundary wall follows the initialization of boundary particles in MPM, meaning that the nodes of the wall condition coincide with the boundary particle positions to avoid errors occurring from data mapping. For the calculation of the contact forces within the DEM partition a Hertzian contact law is considered. Friction and adhesion between the wood clusters itself is set to μ = 0.6 and c = 30 Pa whereas for interaction with the rigid boundary μ = 0.3 and c = 60 Pa is assumed, respectively. For the simulation, a time-step of t = 5e −5 is considered.
The numerical results of the newly proposed partitioned MPM-DEM coupling scheme are presented in Fig. 10 in comparison to the experimental results published by [36]. It can be seen, that our solution can reproduce accurately the experimental result. Similar to the experiment, the granular flow reaches the wooden blocks at t = 0.25s. Due to the impact, the DEM clusters start to move and rotate around the fixed block at the bottom. Comparing the results at time t = 0.30s one can observe, that the two cluster blocks are still in touch over the full width of the block which agrees well with the experiment. This is an enhancement if compared to the numerical results obtained by [37], which is illustrated in Fig. 12, as the adhesion between the DEM blocks can be considered by the implemented Hertzian contact law within the DEM application. In our simulation, the second block firstly touches the ground between t = 0.36s and t = 0.37s with one edge and rotates subsequently until the complete outer edge is in touch with the ground at t = 0.40s which corresponds to a rotation of 90 • in total. Figure 11 shows the rotation angle of the second block in comparison to the experiment and the results from the literature [36,37]. The obtained solution is in good agreement for the entire simulation time.
Block number one, the top block, starts to move and rotate together with the second block and touches the ground for a first time between t = 0.38s and t = 0.39s with its corner. This impact causes its detachment from the ground and its rotation around its axis, which coincides with the picture of the experiment at t = 0.40s and finally comes to rest with a total rotation of 180 • until t = 0.45s. In Fig. 12 the experimental results during the impact are presented in comparison with the available numerical solutions taken from [36] and [37] and our novel approach. In particular, the rotation of the upper block was not captured correctly in previous solutions but is reproduced very well by the novel approach.   [36] and the numerical solutions obtained by [36], [37] and our simulation results during the impact Finally, comparing the run out configuration of the simulation with the experiment, the resting distance of the two moving blocks can be measured. In the experiment [36] the distance between the left boundary of the second block and the left boundary of the box was measured at 34.6cm whereas in our simulation the distance turned out to be 34.5cm. Likewise for the top block, where the distance in the experiment was measured at 40.1cm compared with 39.7cm in the simulation which is in very good agreement. In order to achieve an even higher accuracy of the results, the calculation parameters, especially within the DEM application, have to be calibrated for the respective experiment. Additionally, to prevent the penetration of material points into the DEM particles at the very end of the simulation a finer discretization of the interface as well as a smaller background element size in MPM could be chosen. Apart from that, the numerical solution agrees very well with the experimental results, proving the applicability of the proposed partitioned MPM-DEM coupling scheme for large strain flow events in interaction with discrete obstacles, which is a common case in the numerical investigation of mass-movement hazards.

Conclusion
In contrast to previous solution techniques for the coupling of MPM and DEM, the proposed partitioned coupling scheme offers the possibility to choose the most suitable solver for each involved subsystem including the respective available toolboxes within these solvers. Therefore each physic involved can be solved independently in its preferred reference frame whereas the interaction is transferred to the shared interface.
Assuming mass-movement hazards interacting with discrete obstacles, the continuous flow can be discretized by MPM, which can simulate the large strain flow event also for large scale application cases using a continuum based approach. The discrete obstacles within the flow and the mutual interaction as well as the interaction with the surrounding gravity-driven flow, however, are simulated by DEM, as the contact detection itself and the calculation of the contact forces is advantageous within the DEM application.
During the simulation process, as the MPM and DEM applications are treated as blackbox solvers, boundary conditions are introduced within each partition to ensure the communication of the involved solvers at the shared interface. Within the DEM model, the interface of the surrounding material is represented by a wall condition, which moves according to the material points in MPM. The contact forces between the wall and the DEM particles are mapped to those boundary particles which are discretizing the interface in the MPM partition and are receiving the contact forces as external forces. Due to the applied forces, the material points as well as the boundary particles are moving accordingly, leading to a kinematic update of the shared interface which, in turn, is mapped back to the DEM application to update the wall condition accordingly. The proposed coupling sequence is solved in a weak sense, applying Gauss-Seidel communication pattern.
Several numerical examples have been simulated to asses the quality of the proposed work. The comparison of the obtained results with analytical reference solutions from the literature proves the accuracy of the proposed scheme in both two-and three-dimensional models. Furthermore, the application for elasto-plastic regimes is demonstrated by comparing the results to experimental data from literature. The considered experiment consists of granular material, which flows down due to gravity and impacts three wooden blocks which are placed on top of each other. The numerical solution of this example, where the granular material is modelled by MPM and the wood blocks are discretized by DEM cluster particles, agrees very well with the experimental data for the entire simulation time. Moreover, compared to existing MPM-DEM coupling schemes in the literature, no special treatment such as transforming material points to DEM particles or introducing additional material points at the outline of the DEM domain is required for contact detection as the two involved physics are solved in a partitioned manner.
To conclude, the proposed partitioned MPM-DEM coupling scheme is a powerful method for the simulation of large strain flow events interacting with discrete objects as for example the numerical investigation of mass-movement hazards including huge rocks or other significant obstacles.