A hybrid material-point spheropolygon-element method for solid and granular material interaction
Yupeng Jiang, Minchen Li, Chenfanfu Jiang, Fernando Alonso-marroquin
AA hybrid material-point spheropolygon-element method for solid and granular material interaction
YUPENG JIANG* , MINCHEN LI , CHENFANFU JIANG , FERNANDO ALONSO-MARROQUIN School of Civil Engineering, The University of Sydney, Sydney, Australia 2.
Department of Computer and Information Science, The University of Pennsylvania, Philadelphia, United States
Abstract:
Capturing the interaction between objects that have an extreme difference in Young’s modulus or geometrical scale is a highly challenging topic for numerical simulation. One of the fundamental questions is how to build an accurate multi-scale method with optimal computational efficiency. In this work, we develop a material-point-spheropolygon discrete element method (MPM-SDEM). Our approach fully couples the material point method (MPM) and the spheropolygon discrete element method (SDEM) through the exchange of contact force information. It combines the advantage of MPM for accurately simulating elastoplastic continuum materials and the high efficiency of DEM for calculating the Newtonian dynamics of discrete near-rigid objects. The MPM-SDEM framework is demonstrated with an explicit time integration scheme. Its accuracy and efficiency are further analysed against the analytical and experimental data. Results demonstrate this method could accurately capture the contact force and momentum exchange between materials while maintaining favourable computational stability and efficiency. Our framework exhibits great potential in the analysis of multi-scale, multi-physics phenomena.
Keywords:
Multi-body numerical simulation, Spheropolygon discrete element method, Material point method, Coupling algorithm, Granular mechanics INTRODUCTION
The numerical simulation of the multi-body system is crucial for understanding several key issues in geomechanics, such as the mechanical properties of the complex granular matrix , the interaction between the debris flow and solid structures . It could also benefit the physics-based simulation in computer graphics (CG) to generate photorealistic visual effects for solid-fluid (or granular media) animation . These systems commonly consist of individual bodies with disproportional sizes and various shapes. The fines-ballast granular structure typically exists at the railroad foundation where the ratio between the volume of small (fines) and big particles (ballast) reaches a level of 1:10 . Components in the system could have distinctively different Young’s modulus. For example, the crumb rubbers are artificially added to the railroad or highway foundation for preventing the breakage of surrounding granite particles and further improving the stability . It is also common in CG production to animate intricate multi-body frictional contact between soft and rigid objects. The numerical methods for these issues often require an accurate simulation for the deformation of the individual bodies and contact forces among them. In other ords, both the elastodynamic and the kinetic behaviours of each body need to be properly calculated. The Finite element method (FEM) is commonly applied, and proper treatment of the spatial discretization is essential for each component of the system
11, 12 . However, the accuracy and stability of FEM suffer from the potentially strong distortion of the meshes. Such a problem can be alleviated only by remeshing schemes which, on the other hand, highly compromises the computational efficiency. The material points method (MPM) offers an attractive alternative approach. It discretises the computational domain with meshless particle and therefore avoids the difficulties encountered during large mesh deformation or topology changes. Because the deformation history of the material domain is stored and represented by material points, MPM utilizes a fixed Eulerian background grid that is not distorted during the simulation. MPM and its variants have been successfully applied for the study of continuum granular materials , crack propagation in snow , and computer animation . As a numerical method based on continuum mechanics, either FEM or MPM requires a spatial discretization for every part of a multi-body system. Even for many engineering applications where only the kinematic information of near-rigid bodies is focused, the calculation of continuum partial differential equations still needs to be performed upon them . The high stiffness of the material requires an extremely small time step interval for the stability of the simulation, while the convergence of iterative solvers in the implicit scheme is slow. Meanwhile, both methods use relatively complicated ways to handle the collision and contact force . The computational efficiency could be largely improved if the near-rigid components are properly simplified so that only an accurate description of kinetic information is provided. The discrete element method (DEM) has high efficiency for simulating the Newtonian movements of rigid bodies. It treats near-rigid bodies, or particles, as perfectly rigid and defines an interaction zone that is coated outside each of them. The contact forces are calculated based on the depth, area or volume of the intersections among particle zones. This reasonable simplification ignores the deformation for the particles while accurately simulates the momentum and contact forces through contact laws. No spatial discretization or mesh generation is involved. The coupling between DEM and other numerical methods (FEM, MPM) potentially provides an optimized balance between computational time and accuracy. Several methods have been systematically developed. For example, the DEM-FEM approach is conducted for hierarchical multiscale modelling of granular media and the interaction between the tire tread and granular terrain . The coupling between DEM and Lattice Boltzmann method (LBM) is applied to study the interaction between water and soil particles in hydrological problems . In these studies, discrete elements are simplified to circular/spherical shapes. Results of these simulations well agree with experimental data and analytical solutions. But still, geometrical parameters impose a strong influence on the movement of a near-rigid body. This factor has been further included in the numerical tools like DEM-CFD coupling between fluid and non-spherical particulate system . One of the recent research for coupling of MPM-DEM was conducted by Liu et al . The authors modelled a 2D sand pile collapsing and impacting three rectangular wooden blocks. The granular flow is simulated with MPM. A shrunken-point DEM is used to calculate the movement of blocks. Nine material points are attached at the corners, centre of the edges and the geometrical centre of ach block. The contact between DEM blocks and the MPM flow are detected through the mutual background grid node. Contact forces are calculated at the mutual projection gird based on the momentum information and applied to geometrical nodes on the DEM parts as body forces. The numerical results show a reasonable agreement with the experimental data. This method could further help the damage analysis of buildings under the impact of debris flow. An MPM-DEM hybrid method also developed to exploits the dual strengths of discrete and continuum treatments. However, this method is mainly focused on the improvement of the efficiency of DEM for simulating the granular flow. Discrete elements and material points are replacing each other under certain criterion rather than coexisting and interacting. Liu et al ’ s method unifies the coupling procedures under the computational frame of MPM. It inevitably inherits an issue of MPM for handling the contacts, which is the strong dependency of the background grid system. First, the collision handling in MPM is calculated based on the proportions of momentum that different objects mutually contribute to grids. Contact detection and the calculation of force are fundamentally affected by the resolution and the structure of the grid system. Secondly, the contact algorithm needs to separately calculate the momentum for each object. The simulation will be computationally expensive if it involves a large number of discrete bodies. Meanwhile, the influence of shape is also omitted by simplifying the DEM particles with few material points. For many cases, such simplification is non-trivial and not rigorously discussed. The accuracy of the angular momentum of DEM particles is highly compromised since the forces are only applied at the geometrical nodes. These problems limit the performance of the coupling method and diminish the advantage of using DEM as a rigid body simulator. Therefore, we believe it is necessary to develop an advanced numerical method that could mitigate these issues and provide a more general way for the coupling between DEM and MPM. In this paper, we develop a different coupling method between the DEM and MPM. The rigid body is represented by spheropolygon discrete element method (SDEM) . It was developed for simulating the movement of irregular discrete particles. Comparing with other DEM-based methods, the SDEM could efficiently handle the contact forces among irregular particles and preserve the conservation of mechanical energy. Due to these merits, SDEM has been coupled with the boundary element method (BEM) for simulating sub-particle stress and particle breakage . The interactions between the fluid and irregular rigid bodies are also studied using the LBM-SDEM . In this work, the MPM is used for simulating the deformable part of the multi-body system. Two different constitutive models, linear elasticity and elastoplastic Drucker-Prager model are used to study the elastic-rigid and granular-rigid coupling respectively. The contact between SDEM and MPM are detected with the Euclidean distance instead of the existence of mutual projection grids. Coupling forces are directly calculated with well-established DEM contact model. For rigid particles, forces are applied at the exact contact position. In summary, our method unifies the coupling procedure under the computational frame of DEM. It significantly reduces the coupling dependency to the MPM background grid. The contact detection and force calculation only happen at the boundary of the rigid particles, which are more efficient than that of the pure MPM. The influence of particle shape can be better preserved and no longer needs to be simplified with material points. The paper is organized as follows. The computational methodology for MPM and SDEM is ntroduced in Section 2 and 3 respectively. The coupling method is presented in Section 4. A serial of numerical tests is conducted in Section 5. Results are rigorously discussed with analytical solutions and experimental data as validations of the MPM-SDEM. Potential applications of this method are provided in Section 6. General conclusions are presented in Section 7. The material point method is a hybrid scheme utilizing both Lagrangian particles and Eulerian grids
29, 30 . It follows the governing equations of the continuum mechanics and its discretization is derived from the Galerkin weak form of momentum conservation, similarly to the finite element method. However, unlike the FEM, which discretises the computational area into piecewise subdomains on a mesh, the MPM uses the particle-wise material regions to represent the continuum. Lagrangian variables such as mass, momentum, and position are carried by the material points. The embedding relationship between Eulerian grids and material points is commonly defined by the nodal shape functions with H regularity. At each time step, Lagrangian variables carried by material points need to be first transferred to the corresponding grid nodes. The equation of motion is solved at the grid nodes while volume integrals are approximated through particle quadrature, the velocity is updated accordingly and then interpolated back to material points for their advection and strain updates. Eulerian grids are restored to a standard Cartesian configuration after each time step, only the values and derivatives of the nodal shape functions are constantly recalculated at the beginning of the next time step. Governing equations and discretization for MPM
The Lagrangian kinematic description of a continuum body needs to satisfy a group of partial differential equations (PDEs), including conservation of mass, momentum, and energy. These PDEs are known as governing equations which, combined with the material constitutive model and boundary conditions, determine the behaviour of the material. The conservation of mass is inherently satisfied in MPM since the material points in this study are assigned with constant mass values. The conservation of energy is guaranteed because the simulation assumes an isothermal setting which does not involve the exchange of heat. Therefore, the dynamic state of the material can be obtained by solving the conservation of moment : , ij j i i b u + = , (1) where ρ is the density of the material; u i denotes the displacement, the dots are the notation for the order of time derivative; σ ij is the Cauchy stress tensor, the subscript denotes the components and the deviator of the tensor; b i is the body force term. The PDEs follow Einstein notation. Equation 1 can be solved in the domain Ω through its weak form: * , ( )d 0 i ij j i i u b u + − = , (2) and boundary conditions: ( ) t = t x t , ( ) u = u x u , (3) where the computational domain is denoted as Ω; u * represents the virtual displacement which quals to zero on the boundary section Γ u ; the value of traction t ̄ i is known at the boundary Γ t . Combining the boundary conditions Equation 3 can be further written as: d d d d * * * *i i ij,j i i i i i u u u b u t u + − − . (4) FIGURE 1
Discretization schemes for the material point method As shown in Figure 1, the domain is discretised by the material points (red dots), here also called “Lagrangian Points”. the information of deformation gradient, mass and momentum are carried by these Lagrangian points. Eulerian background grid nodes (blue squares) are defined as background scratchpad. At each time step, the variables are first interpolated to grid nodes using multi-dimensional shape functions. The information is then updated at grid nodes and transferred back to the material points for the next time step. The perspective of generalized interpolation material point method (GIMP) is adopted for the discretization process of the governing equation Eq. (4). Each material point occupies a partition Ω p in the entity Ω: ( )d p p p V = x , (5) where V p is the initial volume and χ p ( x ) is the characteristic function, the subscript p denotes the value on the material point. The mass of the material m p can be written as: ( ) ( )d p p p p m = x x , (6) where ρ p =m p / V p is the density of the material. There are mainly two different forms of χ p ( x ). One is using the total discrete approach (DMPM): ( ) ( ) p p p V = − x x x , (7) where x p is the spatial coordinate of the material point and δ is the Dirac delta function. The mass only exits at the discrete position over the entire computational domain. The other method is to use the GIMP. where the form of χ p ( x ) for each material point to be a continuous function (constant, linear or even a higher-order) over the Ω p . A given physical variable k in the computational domain can be approximated by the value k p carried by the relevant material points and their characteristic functions: ( ) ( ) p pp k k = x x , (8) Equation (4) can be converted from the continuous integration form into a summation of material oints: * * , d = u d p p ip p ip ijp p ip jp pp p uV − * * d d p p p ip ip ip ipp p m b u t uV + + , (9) The behaviours of the continuous material are now defined by the physical variables carried by material points. Equation 8 needs to be solved on the background grid nodes. The relation of virtual displacement between the grid node and material points is written as: * * ( ) p Ip p II N = u x u , * *, , ( ) p j Ip j p II N = u x u . (10) The subscript I is the indexes that denote a value is on the grid node I . Substituting Equation 10 into Equation 9 to eliminate the virtual displacement and the equation of motion can be rewritten as: int ext I I I = + p f f , I u x (11) where p I is the momentum for the grid node; f int I and f ext I represent the internal and external force applied on the grid node respectively: iI ip Ipp p p S = (12) int , iI ijp Ip j pp f S V = − (13) ext d iI p ip Ip i Ip f m b S t N = + (14) The interpolation of variables between the material points and the grid nodes is calculated using the weighting function S Ip ( x ):
1( ) ( ) ( )d p Ip p Ipp
S NV = x x x (15) , ,
1( ) ( ) ( )d p Ip j p p Ip jp
S NV = x x x (16) This interpolation is often referred as P - G (particle to grid) transferring and the reverse process is G - P (grid to particle) transferring. The specific forms of S Ip ( x ) depends on the choice of interpolation method (DMPM or GIMP) and shape function, which is provided in the Appendix. The interpolation function in GIMP has C continuity even if the shape function is only a linear function with C continuity. Generally, GIMP is more stable for spatial discretization, it produces less computational noises when the material points moving from one grid cell to another (cell-crossing noise). However, DMPM with quadratic or cubic B-spline weighting functions also has its advantage for being able to use a noise-free and angular momentum conserving transferring scheme called affine-particle-in-cell method (APIC) . Note that as pointed out by Gao et al ., that GIMP and quadratic-B-Spline-DMPM are equivalent when particle domain is chosen to be a box with its width equal to the grid ell spacing. Time integration scheme
The computational domain is time-variant, which indicates Equation 11 needs to be fulfilled at each time step. The central difference method is applied for the update of the grids’ momentum: n n nI I I t + − = + p p f , (17) where superscript n denotes the sequence of the time step and Δ t is a constant interval of the time increment at each step; f In equals to f I ext + f I in t is the total nodal force. In this paper, we use the update-stress-first (USF) format for the optimal stability and the conservation of energy of the simulation. The algorithm to update the nodal and point variables is introduced in its computational sequences: a ) The nodal mass m I , momentum p I , and velocity v I are calculated through the interpolation of the corresponding material points and their weighting function S Ip : I n np Ipp m m S = , (18) n n nI p p Ipp m S − − = p v , (19) / n n nI I I m − − = v p , (20) b ) The material point’s strain ε ̇ ijpn -1/2 and spin Ψ ̇ ijpn -1/2 rate tensors are obtained with the nodal velocities : n n n n nijp Ip j iI Ip i jII S v S v − − − = + , (21) n n n n nijp Ip j iI Ip i jII S v S v − − − = − , (22) where the value of S Ip is decided according to the method of spatial discretization and the relative position between the material point and grid node. The strain and stress are then calculated as: n n nijp ijp ijp t − − = + , (23) n n nijp ijp ijp t − − = + . (24) The stress rate is determined by: ij ij ik jk jk ik = + + , ij ijkl kl C = , (25) where Jaumann stress rate σ ij ∇ for linear elasticity is adopted to eliminate the influence of pure rotation to the Cauchy stress tensor. The specific form of stress rate σ ij ̇ depends on the constitutive model. c) The internal and external nodal forces are updated using the material point’s stress, the body orce term b p , and the boundary traction t p according to Equation 12 and 13 d) Equation 11 is solved at the grid nodes with the nodal forces and momentum to further update the velocity and position of material points. There are three methods for this specific G-P transferring procedure, particle-in-cell format (PIC) , fluid-Implicit-particle format (FLIP) and the hybrid format. The PIC transferring update the material’s velocity directly using the interpolation functions and corresponding grid values: / n n n np I Ip II S m + + = v p . (26) This transferring method has a strong numerical dissipation; the angular momentum and kinetic energy of material points decrease rapidly with the progress of iterations. This problem makes the PIC format unsuitable for simulating dynamic problems such as the landslide. Because the kinetic energy of granular flows is not accurate enough and hence make the results unreliable. However, this dissipation effect is not entirely undesirable. It efficiently decreases the vibration in quasistatic simulations, which generated by the initial configuration of the objects, and help the system to reach an equilibrium. The FLIP updates the velocity using the nodal force: / n n n n np p I Ip II t S m + − = + v v f . (27) It greatly improves the conservation of angular momentum. But the exact conservation can only be achieved by using the ‘full’ mass matrix instead of lumped mass matrix, which is necessary for the numerical stability but impractical due to its potential singularity . Meanwhile, FLIP format generates the so-called “ringing instability” because of noisy velocity modes in the null space of the transfer operator. It causes the numerical instability when simulating the granular material with strong dynamic behaviours. A better way to preserve the advantages of PIC and FLIP is to use a hybrid transferring format : / (1 )( / ) n n n n n n n np I Ip I p I Ip II I S m t S m + + − = + − + v p v f , (28) where α is the coefficient that ranges between 0 to 1; α= α= α is proportional to the magnitude of the dissipation effect and can be used for controlling the damping for the quasi-static simulation. The instability of FLIP format is also alleviated. However, the hybrid method still has an issue. No rigorous analysis to define the quantitative relation between the dissipation and the value of α . The behaviour of the material is also sensitive to the variation of α and the numerical damping highly dependents on the time step size. Especially for the granular material, the specific value of the coefficient is rather an empirical setting. Such a problem can be properly solved by using affine-particle-in-cell format (APIC) or extended-particle-in-cell format (XPIC) . We adopt APIC in this paper. This innovative method is developed by Jiang and Schroeder . It represents particle velocities as locally affine, which allows APIC to conserve linear and angular momentum across transfers. This transferring format effectively reduces the numerical dissipation; it also does not experience the velocity noise and instability in FLIP. It has been applied for the simulation of both the granular and hyper-elastic material and exhibits superior performance. The APIC still use the PIC format for the velocity and position update of material points at step ( d ). However, APIC applies a different scheme for the P-G transferring procedure at step ( a ), which is written as : / 3( )( ) / 4 n n n n n n Tp Ip I p I pI LN L = − − = ID x x x x I cubicquadartic , (29) ( ( ) ( )) n n n n n n nI p Ip p p p I pp m N − − − = + − p v B D x x , (30) where L is the grid spacing and I is the unit matrix; the additional matrix D p serves as the inertia matrix for affine motion and B p contains the angular momentum information and updates as: ( ) n n n n Tp Ip I I pI S v + + = − B x x (31) In this paper, we use the GIMP with hybrid format for simulating the linear elastic material and DMPM with APIC format for the granular material. The quadratic kernel is adopted as the shape function.
Constitutive models
The physical properties of granular material have been a difficulty for the numerical simulation. Large deformation may happen due to the elastoplastic behaviour of granular material which causes trouble for the mesh-based method. The MPM, therefore, has a unique advantage for simulating both quasistatic deformation and granular flow without the restriction of mesh. The Drucker-Prager plasticity model( D - P ) is employed to simulate the granular material. This model has been widely applied for the engineering application and MPM sand simulation. Although the Jaumann stress rate introduced earlier may not be exceedingly accurate for moderate deformation with deviatoric strains that more than 10 percent (e.g. granular flow) , it can still provide reasonable results and we choose to use it for a direct comparison to the results in the existing literature . FIGURE 2
Yield criteria of Drucker-Prager( D - P ) model The yield criteria in D - P model is shown in Figure 2 and defined as: s f q k = + − , (31) t tm f = − , (32) where f is the yielding surface, the superscripts s and t denote the shear and tensile yielding behaviour respectively; σ t is the tensile strength; τ is the equivalent shear stress and σ m is the spherical stress: J = , /2 ij ij J s s = , (33) /3 m I = , kk I = (34) where J is the second invariant of the deviatoric stress tensor and s ij is the deviatoric stress components; I denotes the first invariant of the stress tensor. The coefficients q ϕ and k ϕ are frictional coefficient and yield stress for shearing behaviour. They are calculated based on the friction angle ϕ and the cohesion term c : q = + ,
39 12 tan ck = + . (35) MPM uses the return mapping method to detect whether the yielding conditions are fulfilled at each time step. The isotropic linear elasticity is used for solid material:
1( )3 ijkl ij kl ik jl il jk ij kl
C K G = + + − , (36) where C is the constitutive tensor; K and G are the bulk and shear module respectively. SPHEROPOLYGON DISCRETE ELEMENT METHOD (SDEM)
The spheropolygon discrete element method (SDEM) is selected as the DEM part of the hybrid algorithm due to its unique advantages. The geometrical irregularity of a particle is properly represented as a Minkowski sum of a polygon with a disk, and the multiple contacts between irregular particles are calculated based on distances between vertices and edges. These features make it computationally efficient . The contact information between MPM and the SDEM can be easily applied to the correct position of SDEM instead of its geometrical modes. The point contact relation is also more realistic for contacts between rigid bodies and granular materials. FIGURE 3
Spheropolygon element obtained by sweeping a disk around a polygon. A spheropolygon is the Minkowski sum of a polygon to represent the irregular shape of near-rigid object and a disk with radius a , which defines an elastic area to calculate the contact forces that generated among particles. Mathematically the Minkowski sum of two sets of points P and Q of a vector space is given by : { | } P Q P, Q + = + x y x y . (37) he geometrical interpretation of this operation is equivalent to the sweeping of a disk around the profile of the polygon while maintaining its original orientation. For example, the DEM particle with a shape of hexagram in Figure 3 is properly approximated by a spheropolygon element with a few boundary lines and a disk sweeping its profile. The hexagram defines the element shape, and the disk is used for contact force calculation. SDEM has been proved as a more effective approach than using a cluster of small particles to approach this shape, especially when the particle shape is more irregular or complicated. The contact force of the spheropolygon is defined through a vertex-edge contact relationship. Let us consider two spheropolygons SP i and SP j with their polygons P i and P j and the radii of the disks a i and a j . Each polygon is defined by its own set of vertices { V } and edges { E }. The overlapping length ξ between each vertex-edge pair ( V, E ) is written as: ( , ) ( , ) i j
V E a a d V E = + − , (38) where d ( V , E ) is the Euclidean distance between the vertex V and the edge E . The brace at the right side of the equation means the non-negative limit of ξ . Therefore, the force vector F applied on particle i by particle j is expressed as: = ( , ) ( , ) i j j i ij ji i j j iV E V E V E V E = − +
F F F F , (39) and the torque τ ij of particle i is: = ( ( , ) ) ( , )+ ( ( , ) ) ( , ) i j j i ij i j i i j j i i j iV E V E V E V E V E V E − − τ p c F p c F , (40) where c i is the centre of mass of particle i and p is the point of contact, which is defined as the middle point of the overlap area between a vertex and an edge:
1( , ) +( ( , ))2 i V E a V E −= − −
X Yp X Y X , (41) where X is the position of the vertex V , and Y is its closest point on the edge E . The movement of the centre of mass r i and the orientation φ i of the particle are governed by the equations of motion: i ijj m = c F , i ijj I = φ τ , (42) where m i and I i are the mass and moment of inertia of the particle; The linear elastic model is used throughout this paper. The force F can be written as: = + n n t t k k F N T , (43) where k n and k t are the normal and tangential stiffness, respectively; N and T are the normal and tangential unit vector that measured at the Edge of the contact; δ n denotes the length of the overlap; δ t is the tangential relative displacement between two particles that accounts for frictional forces. It is limited by the condition k t δ t ≤ μk n δ n , where μ is the friction coefficient. COUPLING of MPM AND SDEM
MPM and SDEM are fully coupled through the contact force information. The crucial part of the coupling is how to properly detect and calculate the contact force between the material points and DEM particles. Liu et al proposed a method which attaches material points to the DEM particles; the square particles are represented with 9 materials at its corners, edges’ middle point and centre of mass. Both the contact detection and contact forces are then calculated using a pure MPM contact algorithm developed by Bardenhangen and Brackbill . This coupling method unifies the contact handling under a pure MPM scheme. The fundamental reason for these problems rises from the simplification of the contact and the dependency of the grid system. The contact force cannot be accurately calculated if only a small amount of material points is involved. However, the advantage of efficiency largely decreases if the material points are densely attached to the DEM particle. Contact handling between MPM and SDEM
In this paper, we proposed a different approach for the computation of contact handling, which unifies the contact detect and force calculation under the scheme of SDEM. The basic idea of this algorithm can be summarized as follows. The Verlet distance, which is the cut-off distance for the potential contact between two discrete elements, is applied to examine the contact between material points and SDEM particles. If a material point is within the Verlet distance of a spheropolygon particle it is treated as a small SDEM disc particle with a certain radius for contact detection. If the material points and spheropolygon is in contact, the contact force between the material point and the spheropolygon is calculated based on DEM contact force models. The calculated force is applied to the material points in a form of extra boundary force term.
FIGURE 4
The contact handling between the spheropolygon and material points (a) Spheropolygon is represented by the green zone and material points by red points, V d is the Verlet radius from the centre of mass o . (b) material points in potential contact with the spheropolygon are assigned with a small radius (red dash circles). As indicated in Figure 4 (a) material points within the Verlet distance V d are transformed into a circular discrete element with its position x p as the centre of mass. Here we name it identified material points (IMP). The interactions among the IMPs are still calculated under the MPM scheme despite the intersections of their radius may happen after the transformation. Figure (b) he magnitude of the contact force f SDEM between an IMP and an SDEM particle can be then calculated based on the modification of Equation 39: ( , ) ( , ) p i p p
E a r d E = + − x x , (44) where r p is the contact radius assigned to an IMP; a i is the sphero-radius of the SDEM particle; d denotes the Euclidean distance between the points x p and the edge E on the SDEM particle. The optimal value of r p is still an open question. Here we provide an estimated interval for r p : p p p l / n r l , (45) where l p is the length of the background grid and n is the average number of the material point in each grid. The value of r p must be smaller than the grid length to control the contact happens within or at the boundary of a cell. Meanwhile, it should be larger than the influential length that a material point representing in the computational domain. If an IMP is in contact with multiple SDEM particles, as shown in Figure 5 (b), the particle-spheropolygon forces are summed together. The single contact of IMP is most likely to happen during the coupling since the radius assigned to the IMP is much smaller than the size of SDEM particles. But in extreme cases, such multiple contact relationship will not affect the stability or the efficiency of the coupling. FIGURE 5
The intersection and contact force calculation between an IMP and SDEM. (a) The calculation of the contact force between one IMP and SDEM particle;(b) the total forces exert on an IMP The contact forces are applied at the corresponding material point as an external boundary force: cont SDEM p k = f f , (46) cont cont ( ) I p Ip pp S = f f x , (47) where the f SDEM is the coupling contact force; f p cont and f I cont is the external force on the IMPs and corresponding nodes respectively. Therefore, the nodal equation of motion (Equation 18) is further modified as: (a) (b) nt ext cont I I I I = + + p f f f . (48) In this way, the influence of the contact force is taken into the calculation of MPM at step ( c ). And through this equation, it could further its influence on the movement of material points in the whole computational area. For the SDEM particle, the contact forces are applied at its centre of mass. The total force and torque exerted on the particle in Equation 39 and 40 are now modified as: SDEM SDEM ( , ) ( , ) ( , ) ( , ) i j j i ij i j j i i p p iV E V E p p
V E V E V E = + + +
F F F f x f x , (49) = ( ( , ) ) ( , )+ ( ( , ) ) ( , ) i j j i ij i j i i j j i i j iV E V E
V E V E V E V E − − τ p c F p c F
SDEM SDEM + ( ( , ) ) ( , )+ ( ( , ) ) ( , ) i p i i p p i i p ip p
V V E E − − p x c f x p x c f x . (50) These equations equal to adding extra contact forces that produced by small circular particles to the large SDEM particle, so that the movement and rotation of the SDEM particles are also affected by the IMP that is in contact with it. Since the algorithm is using an explicit form the critical time step Δ t min needs to be determined for the stability of the coupling. This criterion can be written as: min1min min2 min 2 max n lct mk = (51) where the l min is the minimum length of the background grid and c max is the maximum acoustic velocity of the material in MPM; m min is the minimum mass of the SDEM particle. The coefficient κ and κ are used to further guarantee the stability since the equation for critical time step is obtained based on the linear elasticity, which κ =0.8 and κ =0.1 are commonly used. The contact between material points and the SDEM particle can be fully coupled together through the contact force. There are several advantages to using this scheme. Contact relations are detected with the Euclidian distance between the centre of mass of SDEM particle and the position of a material point. The contact list can be saved as a data structure and reused for the next time step if the displacement of the particles and material points is small. It is more efficient than the condition of mutual grid node where the velocity field on each object is constantly recalculated. Contact forces can be directly calculated with DEM contact models instead of using the grid momentum as an indirect approach. Positions for applying contact forces are not restricted by the grid node or the geometrical nodes of the rigid particles. Irregular shapes are considered, and the angular momentum of the rigid body is better preserved. Coupling method procedures.
The program for MPM-SDEM is developed (MPM-SPOLY) using C++ language for further validations and applications. The coupling algorithm in a time step Δ t can be summarized as follows: i) SDEM : contact detection and force calculation Update the contact list of SDEM particles based on the Verlet distance. Update the contact list of SDEM particles and IMPs. Update Vertices { V } for each particle Update the Vertex-edge contact relations between particles within the contact list Calculate the contact force of SDEM particles and apply the gravity Calculate the f SDEM between SDEM particles and IMPs (ii)
MPM : variables transfer between nodes and points Calculate the nodal mass m I , momentum p I and v I . Calculate the strain and spin rate 𝛆̇ p , Ω ̇ p of material points Calculate the stress σ p of material points (iii) Coupling : Update of material points Calculate the nodal force f I ext , f I int and coupling force f I cont Update the nodal momentum p I Update the position x p and velocity v p of material points (iv) Coupling : Update of the SDEM particles Apply the f SDEM to the SDEM particles Update the velocity, angular momentum Update the position of centre of mass for each SDEM particles VALIDATIONS
A serial of tests is conducted in this section for the validation of MPM-SDEM method. They include three essential parts of the interaction between MPM and SDEM bodies: the conservation of energy, the contact force, and the granular-solid interaction. The first two tests are simulated using linear elasticity as the material property of the MPM; the last test is conducted with Drack-Prager model for the plastic deformation of granular material. The conservation of energy is crucial for the stability of the coupling. Contact forces between the MPM and SDEM bodies need to be correct since the contact handling method is unified under the DEM contact. It is also necessary to examine the plastic behaviour between the MPM and SDEM where the variation of contact relationships is far stronger than that of solid cases. Results are rigorously compared and analysed with analytical solutions. It helps to better understand the advantages and limitations of this method as a general scheme for the coupling between SDEM and MPM.
Conservation of energy
Two tests are designed to investigate the exchanges of momentum and the transferring between the kinetic and gravitational potential energy. The main purpose is to examine whether the energy of the system is increased by the collision between MPM and SDEM. The stability of the simulation can only be preserved if no extra energy is introduced into the system by the coupling algorithm. The other important purpose is to investigate if the coupling could well preserve the conservation f energy and how strong will the coupling affect the system energy. The system still has the numerical dissipation caused by pure DEM and MPM.
Exchange of momentum
The exchange of momentum is conducted by the 2D-elastic collision between MPM and SDEM disc. As illustrated in Figure 7 the MPM disc moving with a constant velocity of 2m/s towards the SDEM body, which is in a static state. The discs have the exact same shape, density, and size. The effect of gravity is eliminated from this test. Material properties and simulation parameters are listed in Table 1.
Table 1
Parameters and material properties for the simulation of conservation of energy
SDEM Parameters MPM (GIMP) Parameters k n Normal stiffness 6.0×10 N/cm d g Grid interval cm k t Tangential stiffness 3.0×10 N/cm r p Coupling radius cm μ Frictional coefficient 0.1 n Number of points 5000 Δ t Time interval 2.0 -4 s MPM material properties V d Verlet distance 0.2 cm ν Poisson’s ratio 0.278 a i Sphero radius 2.0 cm K Bulk modulus / ×10 KPa
SDEM material properties G Shear modulus / ×10 KPa ρ d Density 2.0 g/cm ρ p Density 2.0 g/cm The MPM and SDEM share the same time interval and sequence; the coupling parameters for the collision handling are calculated with the same parameter of normal and tangential stiffness of SDEM. Two groups of material modulus, which are marked bold in Table 1, are used for the MPM disc. The bulk and shear modulus of the hard group is ten times larger than that of the soft case. This comparison is proposed to investigate the influence of the material properties to the conservation of energy. It could further demonstrate the ability for MPM-SDEM to simulate the soft-rigid multibody system.
FIGURE 7
The collision between the MPM and DEM disc at each time slice (Unit: cm/s). (a) the soft material case; (b) the hard material case. (color encodes velocity)
FIGURE 8
The variation of velocity square of DEM and MPM disc with different material properties. It can be observed from Figure 7 (a) that the contact between two bodies is captured by the coupling algorithm and the collision is calculated with the contact force. The velocity of the MPM disc drops rapidly during the collision and transfer its kinetic energy to the DEM disc at t =1.46s and 1.50s. Part of the energy is stored in the MPM body in a form of strain energy because some of the material v x (b) (a) oints at t =1.70s still have a small velocity. This part of the energy is dissipated in the MPM disc since the collision is completed. It dissipated due to the P - G transferring scheme and eventually disappear(e.g. t =2.4s). This phenomenon is consistent with the variation of squared v x in Figure 8, where kinetic energy is largely transferred into the DEM disc, which obtained a v x =1.92m/s. It still shows a clear loss of kinetic energy as these part are not transferred to DEM after the collision. For the hard material case in Figure 8, the conservation of energy is much better than that of the soft case. The collision happened in a short time of duration as it is shown in Figure 7 (b) after t =1.46s. The velocity of material points becomes zero after the collision. The SDEM disc obtained a v x =1.992m/s. The kinetic energy is much closer to the analytical solution of the perfect elastic collision. Transferring between the gravitational potential and kinetic energy
The test for the potential-kinetic energy transferring is conducted by dropping an MPM elastic disc at the spheropolygon boundary. The spheropolygon element is set as an unmoveable elastic boundary. It can be seen from Figure 9 that the MPM disc hits the boundary and bounces back at the spheropolygon boundary. The contact is detected by the algorithm and the collision are properly calculated. The simulation parameters and material properties are the same as Table 1; the gravitational acceleration is -10.0 cm/s FIGURE 9
The soft elastic disc dropping at and bouncing back by the spheropolygon boundary The variations of the vertical velocity square are presented in Figure 10. The hard disc is shifted with +0.5s on the time axis for the comparison. Figure 10 indicates that the soft discs still have a clear loss of energy after the collision. The hard disc case, similarly, has a better performance for the conservation of energy; less energy is dissipated during the collision; its velocity almost reaches the theoretical limit. It can also be observed from Figure 11 that soft disc experienced a longer time to complete the transferring between kinetic and potential energy. For both soft and hard disc, no extra energy is generated by the collision since the squared velocity after the collision is lower than the analytical limit. It indicates the coupling algorithm can maintain its stability. This results further supports the fact that the material property has a certain influence on the conservation of energy. The low value of elastic modulus could generate a strong loss of energy after the coupled collision. It may cause a problem for the simulation that requires a high velocity (e.g. impact engineering) but could still be well-applied for the simulation of quasi-static or low-velocity issues.
FIGURE 10
The change of velocity squared of the MPM disc; the hard disc is shifted with a 0.5s on the time axis for the comparison.
FIGURE 11
The time gap t and t of collision obtained with soft and hard MPM disc case. Coupling force
The coupling of the MPM and SDEM is conducted through the contact force that calculated based on the contact model of DEM. It distinguishes the coupling method developed in this paper from the existing method, which can be used for the quasi-static and dynamic cases. The value of normal and frictional coupling force between the MPM and SDEM are examined in this section. Especially, the contact force’s dependence on the grid size is investigated here with four different values of grid interval. Simulation parameters and material properties are given in Table 2; integer i denotes the variation of the grid size, which ranges from 0 to 3 and marked in bold. Table 2
Parameters and material properties for the simulation of MPM-SDEM contact force
SDEM Parameters MPM (GIMP-PIC) Parameters k n Normal stiffness 6.0×10 N/cm d g Grid interval 0.35+0.05× i cm k t Tangential stiffness 3.0×10 N/cm r p Coupling radius cm Frictional coefficient 0.1 n Number of points 2601 Δ t Time interval 2.0 -4 s MPM material properties V d Verlet distance 0.2 cm ν Poisson’s ratio 0.2558 a i Sphero radius 0.5 cm K Bulk modulus 6.0×10 KPa General Parameter G Shear modulus 3.5×10 KPa g Gravitational acceleration 100.0 cm/s ρ p Density 2.5 g/cm Normal force
The normal contact force is examined by placing an MPM elastic square on the spheropolygon boundary. Reaction forces are generated due to the gravitational force and, based on the coupling algorithm, applied on each material point at the bottom of the square. Figure 12 indicates that the contact relations between material points and the boundary are properly captured; the contact force is equally applied at the points that are in contact with the spheropolygon boundary. There are in total 51 material points near the boundary and each one has a 0.4902N normal force as the reaction force for the gravity.
FIGURE 12
Normal forces generated between the elastic MPM block and spheropolygon boundary (Unit: N); d g =0.35 and t =3.5s.The normal forces are generated at the bottom of the MPM elastic square where the material points are in contact with the spheropolygon boundary; the magnitude of the force generated each point are all equal to 0.4902N (red value in colour bar); the other material points have zero force value (blue value in colour bar) since they are not in contact with the boundary. The total reaction force for different grid size is shown in Figure 13. It indicates that the simulated value of normal force has a strong vibration at the beginning of the simulation. Such vibration is caused by the initial zero-overlapping configuration between the boundary and the elastic block. The normal force gradually converges to the analytical solution. The progress indicates that the static equilibrium is searched by the contact algorithm and eventually reached. The kinetic energy is dissipated due to the viscoelastic contact model of DEM and the numerical damping within the point-grid transferring scheme of MPM. The energy dissipation here is a positive factor, which helps the system reaches a static state and provides the correct force information. This tendency can be observed for all four cases with different grid sizes. The larger size of the grid only improves the intensity of the energy dissipation but does not change the convergence of the contact force; the (a) (b) orrectness of the contact force is independent of the size of the background grid. The MPM and SDEM can be properly coupled through the DEM contact scheme. FIGURE 13
The total normal force for different background grid size.
Frictional force
The same model is used for simulating the frictional force. A constant velocity of V x =2.0m/s is applied to the same MPM square. It is sliding towards the positive direction of the x -axis and generates a friction force. Its friction coefficient between the MPM and DEM boundary are set as μ= μF n =- FIGURE 14
Frictional forces generated between the elastic MPM block and spheropolygon boundary (Unit: N); d g =0.35 and t =3.5s. The frictional forces are generated at the bottom of the MPM elastic square where the material points are in contact with the spheropolygon boundary while sliding; the magnitude of the frictional force on each point are all equal (blue value in colour bar) since their relative position to the boundary is the same; the other material points have zero frictional force (red value in colour bar) since they are not in contact with the boundary. The variation of total frictional force is shown in Figure 15. The test for the effect of grid size has been clarified in the normal force test and therefore not repeated here. The grid size d g =0.50 is used (a) (b) o quickly reduce the vibrations generated by initial configuration. Similarly, the friction force also converges to the analytical value with the progress of the simulation. The vibration energy is dissipated and the MPM square reaches a quasistatic state; both normal and friction force are correctly calculated. FIGURE 15
The total frictional force for d g =0.5 Granular flow
The coupling performance between the SDEM and the MPM granular material are tested with the silo flow model. The Drucker-Prager model is used to simulate the plastic behaviour of the non-cohesive dry sand which represented with MPM. The SDEM are used as the silo and container of the granular flow. Five different diameters of the silo neck are used to compare the results with the 2D Beverloo law. The main reason for conducting these tests is because the contact relationship for the plastic MPM-SDEM case is far more changeable than that of elastic MPM-SDEM cases. Contacts are constantly generated and deleted in the granular flow. Therefore, it is necessary to test whether coupling simulation is reliable. Material points, although representing a continuous area, are used to approximate the assemble of discrete bodies (i.e. sand or rock pile). Three different number of material points are used to investigate its influence on mass transportation. The simulation parameters and material properties are given in Table 3.
Table 3
Parameters and material properties for the simulation of MPM-SDEM granular flow
SDEM Parameters MPM (APIC) Parameters k n Normal stiffness 6.0×10 N/cm d g Grid interval 0.35 cm k t Tangential stiffness 3.0×10 N/cm r p Coupling radius cm μ Frictional coefficient 0.2 n Number of points 8/9/10×10 Δ t Time interval 2.0 -4 s MPM material properties V d Verlet distance 0.2 cm ν Poisson’s ratio 0.2358 r d Sphero radius 0.5 cm K Bulk modulus 5.0×10 KPa General Parameter G Shear modulus 3.2×10 KPa g Gravitational acceleration 10.0 cm/s ρ m Density 1.5 g/cm Drucker-Prager model σ t Tensile strength 0.0 MPa ϕ Friction angle 35.0 degree ψ Dilation angle 25.0 degree
FIGURE 14
The v for the dry sand silo flow at D =0.04m (Unit: cm /s )(a) t =0.0s (b) t =1.0s (c) t =2.0s (d) t =3.0s As illustrated in Figure 14, the MPM dry sand flows through the bottleneck of the silo and reaches a steady flowing state; the sand drops at the bottom and held by the SDEM container. It indicates that the coupling algorithm remains stable for the granular material. The contact can be effectively detected and properly handled. FIGURE 15
The mass transfer rate for different silo neck diameter and number of material points; C =0.58 and k c =2.2 are the empirical coefficient; D is the diameter of the silo neck. The mass transfer rates for different sizes of silo neck are shown in Figure 15. Simulations with different numbers of material points are performed for each diameter. Results indicate that the simulated mass transfer rates agree well with the analytical solution of the 2D Beverloo law . The rate increases with the increase of the diameter in a 3/2 power law relation. The number of material points does not have a strong influence on the results since MPM is a continuous numerical method. Although there is still a deficiency for using the Drucker-Prager model as indicated earlier, the coupling method is reliable for simulating granular-solid interaction. Discussion
The studies in this section prove that the coupling between MPM and SDEM can be properly calculated under the contact scheme of DEM. Results agree well with the analytical solutions for tests of conservation of energy, contact force, and granular flow. The conservation of energy is affected by the material property but properly maintained in general. No extra energy is generated from the collision between material points and SDEM contact layer, which indicates the stability of the method. Contact forces, as the communicating value for the coupling, are correctly provided. Such correctness guarantees the proper interaction between a continuous (MPM) and a discontinuous (SDEM) method. Tests of granular flow indicate that the coupling could effectively detect and handle the highly changeable contact relations between the granular media and SDEM particles. Therefore, MPM-SDEM coupling method can be applied for many numerical simulations, which have a large difference in size or Young’s modulus. APPLICATIONS OF MPM-SDEM
Two applications of MPM-SDEM coupling method are provided in this section to further demonstrate its ability to simulating the problems that regarded as challenging for traditional methods. In the first application, we simulated the motion of wooden blocks under the impact of granular materials. The simulation is conducted with the same parameters and configurations used y Liu et al for the coupling of MPM and DEM. It is a meaningful comparison to study the difference between two coupling methods. The second application is a quasi-static compression of a fine-particle granular matrix. It is a typical structure for Geomechanics that has exhibited complication mechanical features and poorly studied due to the limitation of numerical and experimental tools. The motion of wooden blocks under the impact of the granular flow is simulated in two dimensions. This test was originally conducted by Liu et al for the validation of a coupled MPM-DEM method. In their method, each DEM block is represented by 9 visual material points to couple its interaction with MPM granular flow based on the exchange of momentum. The configurations of the test are given in Figure 16. The granular flow will be released and impact three piled-up wooden blocks. For the granular material, the gravitational potential is transferred to the kinetic energy and eventually exerts on the blocks through the point-grid projection on the mutual background grids. Two upper blocks will be pushed to the right side and obtain angular momentum; Block.3 is glued to the ground and thus unmovable. The variation of the inclining angle of block No.2 is recorded and compared with experimental data. FIGURE 16
The initial configuration of the granular flow and SDEM blocks In this section, this test is reconducted using our MPM-SDEM coupling scheme. The results are compared and rigorously analysed with both numerical and experimental data from the existing. It will further illustrate the advantages of using the contact force as the coupling connection. The material properties and simulation parameters are given in Table.3 with a few modifications for the simulation: g =10m/s, ψ =0°, n =8000 and 𝜙 =22.0°. The movements of granular flow and blocks for the simulation are shown in Figure 17. The time for each screenshot is taken at the exact same time of Liu et al’s study to provide a clear visual comparison. They indicate that both the collapse of the granular pile and the motion of the blocks agree well with the existing numerical and the experimental data at each compared time step. The sand flow reaches the pile of blocks at t =2.5s; the blocks start to move and rotate due to the impact. he block No.2 touches the ground at t =4.0s and block No.1 touches the ground after t =4.5s. FIGURE 17
Comparison among the Liu et al’s experiment(left) and numerical results(middle), and MPM-SDEM (right) simulation of granular flow impacting blocks at different time steps. The rotation of block No.2 is recorded and compared with the numerical and experimental data in Figure 18. It can be observed that comparing with the existing method inclination angles of MPM-SDEM are closer to the experimental results at each measured time step, except t =0.3s. Such better performance is particularly obvious at t =0.25 and t =0.40s. The differences between MPM-SDEM and experimental data are smaller than 4.0°, where the different almost reaches 20.0 ° for the existing method. Furthermore, the variation curve of the experiment is smoother than that of Liu’s method. It indicates that the actual rotation and changes of angular momentum of block No.2 are relatively continuous. The result of the existing coupling method exhibits jumps of the angular velocity. This discontinuity is particularly strong between t =0.2s to 0.25s. As discussed earlier, this problem is caused by two aspects. First, the DEM blocks are insufficiently represented with only 9 material points. There is a strong loss of accuracy when calculating the momentum transferred from granular flow to the highly simplified DEM block. Secondly, the coupling method has a strong dependency on the background grid. The angular momentum will be highly compromised if the gird size is too large; however, if the grid size is too small the material point at the centre of mass can be ignored and receives zero momentum from the granular flow. The grid size has a strong influence on the simulation, yet the optimal value is very hard to be determined; neither increase nor decrease the grid size can guarantee the increase of accuracy. The improvement can be made by increasing the number of material points to represent the DEM block. But that strategy makes a large compromise of computational efficiency; both detection and calculation of the contact become far more expensive. Therefore, this method is only reliable for certain cases. The variation curve provided by MPM-SDEM shows a continuously changing pattern like the experimental data. No sudden jumps of angular velocity for the SDEM block. Because the coupling is conducted through the contact force. Neither contact detection nor force calculation has a strong dependency on the grid. The angular momentum is better preserved during the impact process since the transferring of the momentum is not limited by the grid. Material points of the granular flow can transfer its kinetic energy to the edges of blocks at any position if their Euclidian distance is recognized as in contact. Contact detection and force can be obtained using solely the Euclidian distance. MPM-SDEM shows a superior performance comparing with the existing method and hould be reliable for a variety of application due to its unique advantages. FIGURE 18
The temporal variation of block No.2’s inclining angle β. An important application for the MPM-SDEM coupling method is the simulation for granular matrix. It is consists of particles in very different sizes and thus produces different physical properties as the assembling of them. The size ratio between large particle and fines could reach 10 level. It is highly computational expensive to simulate this structure with pure DEM method. Pure continuous methods also have difficulty to properly simulate the interaction between large particles and fines. Studies
1, 44 have been done with DEM in 2D and certain progress has been made based on the results. However, both particles and the fines are poorly represented with a limited number and therefore compromises the accuracy of the results.
The MPM-SDEM method provides a better way to simulate the micro-mechanism of the granular matrix. It is a scale-crossing method that could represent the particles in an extreme size ratio. Large particles are represented by SDEM, the irregular shape can be properly considered. Fines, which exists in the voids between the particles, are simulated with MPM. The mechanical properties of the fines are governed by the D - P model. Its effects no longer need to be represented with a large amount of the tiny particles while the discrete interaction of large particles can still be preserved. An example of the uniaxial compression test is conducted to investigate the effect of fine percentage in the granular matrix. Large particles are generated with irregular pentagons. The parameters are the same as Table 4 except the density of SDEM particles is changed to 2.2 g/cm . FIGURE
19 The uniaxial compression of granular matrix consist of fines and particles. The container size of 61×61 cm; the compression bar is moving with a constant velocity of v y =-0.5cm/s. (a) The initial configuration of the granular matrix (b) granular matrix at t =22.0s The basic structure of the granular matrix sample is shown in Figure 19 SDEM particles are the skeleton while MPM fines exist in the voids among them. Both of the components bear the compressive load from the top. At the early stage, the compression mainly exerts on the SDEM particles since the void space is not fully saturated. SDEM particles will start to slowly move and rotate. It can be regarded as a rearrangement process for the skeleton part of the granular matrix. The void space is shrinking due to this rearrangement. MPM fines start to carry the compressive load by interacting with SDEM particles if the void space is small enough. This phenomenon is one of the reasons that the granular matrix exhibits a relatively complex mechanical behaviour.
FIGURE
20 The force-settlement curve for different percentage of fine content. The compression force-settlement curves are shown in Figure 20. It can be observed from the results that the fine percentage has a major impact on the compressibility of the granular matrix, which would further affect the result for measuring Young’s modulus. The basic skeleton structure of the SDEM is preserved because the void space is not saturated by fines and the skeleton structure thus should not be strongly affected by a small change of fine content. It indicates that the fine content (a) (b) ould effectively change the mechanical property; a higher fine percentage largely reduces the compressibility of the granular matrix and, therefore, exhibits a higher Young’s modulus at the macroscale. The test provided a good demonstration of the potential application of MPM-SDEM method. More rigorous numerical tests with comprehensive consideration and higher accuracy will be conducted in the future study. CONCLUSIONS
An advanced MPM-SDEM coupling method is developed for simulating the multi-body system with a large ratio of size and material modulus. Especially the interaction between solid and granular materials. This continuous-discontinuous algorithm has advantages of the MPM for dealing with large deformation and the efficiency of SDEM for handling the movement of the irregular near-grid body. This method is conducted through the contact force between SDEM particles and material points. Both contact detection and the force calculation are unified under the contact scheme of the SDEM method. A material point is identified as a circular discrete element (IMP) if its Elucidian distance to the centre of mass of an SDEM particle is smaller than the Verlet distance. The normal and tangential forces of the contact can then be calculated based on the intersection and relative displacement between the SDEM particle and the IMP. For MPM simulation, this contact force will be applied to the IMP as an external boundary force term. The movement of material points is influenced and constrained by the contact forces, which is also added to the total force vector on the corresponding SDEM particle to determine its velocity, rotation, and position. A serial of validation tests is conducted with two different constitutive models, linear elasticity and Drucker-Prager model. The conservation of energy, contact force, mass transfer rate in granular flow is investigated and analysed with analytical solutions. Results indicate that the conservation of energy is well preserved by the coupling method through the contact force. However, a certain amount of energy will lose during the collision if the material modulus of MPM is extremely low. Both normal and friction force between the MPM and SDEM particle are correctly calculated after the vibration generated by the initial configuration. Such vibration energy can be effectively dissipated with a larger size of the background grid and transferring format. Granular flow simulated by MPM also provides a good agreement with the Beverloo law for the mass transfer rate. Two representative applications are presented. The motion of wooden blocks under the impact of the granular flow is conducted and compared with an existing MPM-DEM coupling method and experimental data. The variation of inclination angle provided by MPM-SDEM is more continuous and shows a better agreement with the experimental data. It overcomes the problems of the grid dependency in the existing MPM-DEM algorithm. Uniaxial compression tests for a granular matrix are performed. It further demonstrates that MPM-SDEM coupling method provides a scale-crossing solution for simulating the complex mechanical behaviour of the granular matrix. Both continuous behaviour of fines and discrete features of large particles can be well preserved in the simulation with optimal efficiency. The results indicate that the fine percentage in a granular matrix has a major influence on its macroscale mechanical properties. In conclusion, the MPM-SDEM has its unique advantage for simulating the multi-body interaction with highly different size and material modulus. We anticipate it would make a good contribution to the study of geomechanics and CG production.
EFERENCE Andrade JE, Avila CF, Hall SA, Lenoir N, Viggiani G. Multiscale modelling and characterization of granular matter from grain kinematics to continuum mechanics.
J Mech Phys Solid . 2011; 59: 237-250. 2.
Salot C, Gotteland P, Villard P. Influence of relative density on granular materials behavior: DEM simulations of triaxial tests.
Granul Matter . 2009; 11: 221-236. 3.
Kruyt NP. Micromechanical study of fabric evolution in quasi-static deformation of granular materials.
Mech Mater . 2012; 44: 120-129. 4.
Utili S, Zhao T, Houlsby GT. 3D DEM investigation of granular column collapse: Evaluation of debris motion and its destructive power.
Engrg Geol . 2015; 186:3-16. 5.
Dai Z, Huang Y, Cheng H, Xu Q. SPH model for fluid-structure interaction and its application to debris flow impact estimation.
Landslides . 2017; 14: 917-928. 6.
Teufelsauer H, Wang Y, Pudasaini SP, Borja RI, Wu W. DEM simulation of impact force exerted by granular flow on rigid structures.
Acta Geotech . 2011; 6: 119-133. 7.
Gao M, Tampubolon AP, Jiang C, Sifakis E. An adaptive generalized interpolation material point method for simulating elastoplastic materials.
ACM Trans Graph . 2017; 36(6): 233:1-233:12. 8.
Tampubolon AP, Gast T, Klár G, Fu C, Teran J, Jiang C, Museth K. Multi-species simulation of porous sand and water mixtures.
ACM Trans Graph . 2017; 36(4): 105:2-105:11. 9.
Huang H, Tutumluer E, Dombrow W. Laboratory characterization of fouled railroad ballast Behavior.
Transp Res Rec . 2009; 2117: 93-101. 10.
Sol-Sánchez M, Thom NH, Moreno-Navarro F, Rubio-Gámez MC, Airey, GD. A study into the use of crumb rubber in railway ballast.
Constr Build Mater . 2015; 75: 19-24. 11.
Kwan SHJ, Sze HYE, Lam C. Finite element analysis for rockfall and debris flow mitigation works.
Can Geotech J . 2019; 56: 1225-1250. 12.
Imseeh WH, Alshibli KA. 3D finite element modelling of force transmission and particle fracture of sand.
Comput Geotech . 2018; 94:184-195. 13.
Lee NS, Bathe KJ. Error indicators and adaptive remeshing in large deformation finite element analysis.
Finite elem anal des . 1994; 16: 99-139. 14.
Bardenhagen SG, Brackbill JU, Sulsky D. The material-point method for granular materials.
Comput Methods Appl Mech Engrg . 2000; 187: 529-541. 15.
Gaume J, Gast T, Teran J, Herwijnen AV, Jiang C. Dynamic anticrack propagation in snow.
Nat Commun . 2018; 9: 3047. 16.
Wolper J, Fang Y, Li M, Lu J, Gao M, Jiang C. CD-MPM: Continuum damage material point methods for dynamic fracture animation.
ACM Trans Graph . 2019; 119: 1-15. 17.
Bardenhagen SG. Guilkey JE, Roessig KM, Brackbill JU, Witzel WM, Foster JC. An improved contact algorithm for the material point method and application to stress propagation in granular material.
Comput model Engrg Sci . 2001; 2(4): 509-522. 18.
Hu W, Chen Z. A multi-mesh MPM for simulating the meshing process of spur gears.
Comput Struct . 2003; 81: 1991-2002. 19.
Zhang HW, Wang KP, Chen Z. Material point method for dynamic analysis of saturated porous media under external contact/impact of solid bodies.
Comput Methods Appl Mech Engrg . 2009; 198: 1456-1472. 0.
Cundall PA, Strack ODL. A discrete numerical model for granular assemblies.
Geotechnique.
Michael M, Vogel F, Peters B. DEM-FEM coupling simulations of the interactions between a tire tread and granular terrain.
Comput Methods Appl Mech Engrg . 2015; 289:227-248. 22.
Feng YT, Han K, Owen DRJ. Coupled lattice Boltzmann method and discrete element modelling of particle transport in turbulent fluid flows: computational issues.
Int J Numer Methods Eng . 2007; 72: 1111-1134. 23.
Zhong W, Yu A, Liu X, Tong Z, Zhang H. DEM/CFD-DEM Modelling of non-spherical particulate systems: theoretical developments and applications.
Powder Tech . 2016; 302: 108-152. 24.
Liu C, Sun Q, Zhou GGD. Coupling of material point method and discrete element method for granular flows impacting simulations.
Int J Numer Methods Eng . 2018;115: 172-188. 25.
Yue Y, Smith B, Chen PY, Chantharayukhonthorn M, Kamrin K, Grinspun E. Hybrid grains: adaptive coupling of discrete and continuum simulations of granular media.
ACM Trans Graph.
Alonso-marroquín F, Wang Y. An efficient algorithm for granular dynamics simulations with complex-shaped objects.
Granul Matter . 2009; 11: 317-329. 27.
Jiang Y, Herrmann HJ, Alonso-Marroquin F. A boundary-spheropolygon element method for modelling sub-particle stress and particle breakage.
Comput Geotech . 2019; 103087. 28.
Galindo-Torres SA. A coupled discrete element lattice Boltzmann method for the simulation of fluid-solid interaction with particles of general shapes.
Comput Methods Appli Mech Engrg . 2013; 265: 107-119. 29.
Sulsky D, Zhou SJ, Schreyer HL. Application of a particle-in-cell method to solid mechanics.
Comput Phys Commun . 1995; 87: 236-252. 30.
Jiang C, Schroeder C, Teran J, Stomakhin A, Selle A. The material point method for simulating continuum materials.2016:
ACM SIGGRAPH 2016 Courses : 1–52. 31.
Gan Y, Sun Z, Chen Z, Zhang X, Liu Y. Enhancement of the material point method using B-spline basis functions.
Int J Numer Methods Eng . 2018; 113: 411-431. 32.
Bardenhagen SG, Kober EM. The generalized interpolation material point method.
Comput Model Engrg Sci . 2004;5(6): 477-496. 33.
Jiang C, Schroeder C, Teran J. An angular momentum conserving affine-particle-in-cell method.
J Comput Phys . 2017; 338:137-164. 34.
Bardenhagen SG. Energy Conservation Error in the Material Point Method.
J Comp Phys . 2002;180: 383–403. 35.
Evans MW, Harlow F. The particle-in-cell method for hydrodynamic calculation.
Los Alamos Sci Lab report . 1957. 36.
Brackbill JU, Ruppel HM. FLIP: A method for adaptively zoned, Particle-in-cell calculations of fluid flows in two dimensions.
J comput Phys . 1986; 65: 314-343. 37.
Love E, Sulsky DL. An unconditionally stable energy-momentum consistent implementation of the material-point method.
Comput Methods Appl Mech Engrg . 2006; 195: 3903-3925. 38.
Zhu Y, Bridson R. Animation sand as fluid.
ACM Trans Graph . 2005; 24(3): 965-972. 39.
Hammerquist CC, Nairn JA. Anew method for material point method particle updates that reduces noise and enhances stability.
Comput Methods Appl Mech Engrg . 2017; 318: 724-738. 40.
Klár G, Gast T, Pradhana A, Fu C, Schroeder C, Jiang C, Teran J. Drucker-Prager elasto-lasticity for sand animation.
ACM Trans Graph . 2016; 35(4): 103:1-103:12. 41.
Molenkamp F. Limits to the Jaumann stress rate.
Int J Numer Anal Methods Geomech . 1986; 10: 151-176. 42.
Courant R, Friedrichs K, Lewy H. On the partial difference equations of mathematical physics.
IBM J.
Mankoc C, Janda A, Arévalo R, Pastor JM, Zuriguel I, Garcimartín A, Maza D. The flow rate of granular materials through an orifice.
Granul Matter . 2007; 9(6): 407-414. 44.
Fu P, Dafalias YF. Study of anisotropic shear strength of granular materials using DEM simulation.
Int J Numer Anal Methods Geomech . 2011; 35: 1098-1126.
APPENDIX
The three-dimensional weighting functions for the GIMP is written as: ( ) ( ) ( )
Ip xIp p yIp p zIp p
S S S S = x x x , (A1) where S iIp , which is calculated based on the constant characteristic function and linear shape function, is defined as:
22 2 2 p ip iI pip iIiIp ip iI p pip iIp ip iI p
L l x x Llx x LS x x l Llx x LL l x x Ll + + − + −= − − + + − + − − ( )( ) ip iI pp ip iI pp ip iI pp ip iI pp ip iI pp ip iI pip iI p x x L lL l x x L lL l x x ll x x ll x x L lL l x x L lx x L l − − +− + − − +− + − −− − − −− − +− + (A2) where L is the spacing of the background grid and l p is the half size of the square influential area defined by the characteristic function χ p . This weighting function is used for all the simulations of solid material. Since the characteristic function is Dirac delta function χ p ( x )= δ ( x p ), the three-dimensional weighting function for the APIC format is also the shape function which is using the quadratic kernel: ( ) ( ) ( ) Ip xIp p yIp p zIp p
N N N N = x x x , (A3) where N iIp is written as: (1 / 2)( /L) ( / ) (2 / 3),(1 / 6)(2 / ) ,0, ip Ip ip IpiIp ip Ip x x x x LN x x L − − − += − − ip Ipip Ipip Ip x x Lx x Lx x L − − − − − −