# A partitioned material point method and discrete element method coupling scheme – Advanced Modeling and Simulation in Engineering Sciences

#### ByVeronika Singer, Klaus B. Sautter, Antonia Larese, Roland Wüchner and Kai-Uwe Bletzinger

Aug 16, 2022

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 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,text {cm}) and a square cross-section of (text {A}=1 times 1 ,text {cm}^2) whereas the impacting sphere has a radius of (R=1.0,text {cm}) and an imposed velocity of 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,text {cm}) meshed by (62 times 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 times 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 times 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,text {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 between DEM and MPM is considered. For both examples, a time step of (Delta text {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,text {cm}) and the radius of the sphere being (R = 2.0, text {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 times 13) quadrilateral elements considering a total length of (30.96,text {cm}) for the grid. The material points are initialized at the Gauss point positions of a structured body mesh with (24 times 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.

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 times 20 times 20 , text {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,text {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,text {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,text {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 (delta ) 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 (delta ), 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 times 1.8 times 19.8,text {cm}^3) consists of (8 times 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 (mu = 0.6) and (c=30, text {Pa}) whereas for interaction with the rigid boundary (mu = 0.3) and (c=60,text {Pa}) is assumed, respectively. For the simulation, a time-step of (Delta t= 5text {e}^{-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.25text {s}). 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.30text {s}) 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.36text {s}) and (t=0.37text {s}) with one edge and rotates subsequently until the complete outer edge is in touch with the ground at (t=0.40text {s}) which corresponds to a rotation of (90{^circ }) 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.38text {s}) and (t=0.39text {s}) 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.40text {s}) and finally comes to rest with a total rotation of 180(^{circ }) until (t=0.45text {s}). 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.

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.6 text {cm}) whereas in our simulation the distance turned out to be (34.5 text {cm}). Likewise for the top block, where the distance in the experiment was measured at (40.1 text {cm}) compared with (39.7 text {cm}) 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.