MatDEM-fast matrix computing of the discrete element method☆

Chun Liu,Hui Liu,Hongyong Zhang

School of Earth Sciences and Engineering,Nanjing University,Nanjing,210023,China

Keywords:Discrete element method High-performance MatDEM Matrix computing

ABSTRACT Discrete element method can effectively simulate the discontinuity,inhomogeneity and large deformation and failure of rock and soil.Based on the innovative matrix computing of the discrete element method,the highperformance discrete element software MatDEM may handle millions of elements in one computer,and enables the discrete element simulation at the engineering scale.It supports heat calculation,multi-field and fluidsolid coupling numerical simulations.Furthermore,the software integrates pre-processing,solver,postprocessing,and powerful secondary development,allowing recompiling new discrete element software.The basic principles of the DEM,the implement and development of the MatDEM software,and its applications are introduced in this paper.The software and sample source code are available online (http://matdem.com).

1.Introduction

The discrete element method (DEM) was first proposed by Cundall and strack,(1979) to study the motion and interaction of elements.In recent years,with the development of discrete element theory and the advancement of computer technology,the DEM has been widely used in various fields (Xu et al.,2003),including the simulation of grain accumulation in agriculture (Jiang et al.,2014),crushing in mining,metallurgy and grain mixing in industrial production,researching on the mechanical properties of grain accumulation in solid mechanics (Goldenberg and Goldhirsch,2005),and the formation and evolution of various geological structures (Antonellini and Pollard,1995;Hardy and Finch,2006;Hardy and Finch,2010;Liu et al.,2015).

The DEM can effectively simulate geological disasters such as landslides or collapses,and provide a reference for predicting the possibility of disasters and evaluating the possible damage (Jiao et al.,2000).The method has been used to study the law of seismic wave propagation and the impact of earthquakes (Feng et al.,2012;Shi et al.,2013).By applying specific triaxial stress to rock,DEM accurately simulates the deformation and damage that occurs inside the rock under considerable pressure (Li et al.,2018).At both macro and micro scales,the DEM is used to simulate test processes employed in rock and soil mechanics.It can therefore quantitatively analyze microstructure and micromechanical mechanisms under?the macroscopic deformation and failure of rock and soil mass (Li et al.,2003;Li and Wang,2004;Jiang et al.,2012,2014).The DEM provides an effective approach for exploring microscopic mechanisms under?the macroscopic mechanical properties of rock and soil mass.

The major commercial software for discrete elemental numerical analysis worldwide includes PFC (Zhou et al.,2016) and EDEM (Wang et al.,2016).Open source software includes Yade,Esys-Element,and LIGGGHTS (Kloss and Goniva,2011).PFC is a calculation software developed by the Itasca Company of the United States and is mainly used in the fields of mining,geotechnical,and earth sciences.EDEM is the flagship product of DEM Solutions in the United Kingdom.It is mainly used for the simulation analysis of grain processing and production operations in industry.Yade is the most widely used community-driven open source software and is written in C++and Python.After more than a decade of development,many excellent discrete element software programs have been designed in China,including 2D-Block,GDEM,SDEM and DICE2D.2D-Block is a two-dimensional discrete element software program (Wang,2000) developed based on Visual C++6.0.GDEM is a software jointly developed by the Joint Laboratory of Discontinuous Media Mechanics and Engineering Disasters,Chinese Academy of Sciences and GDEM Technology,Beijing Co.,Ltd.Based on the DEM of continuum mechanics,GDEM realizes the entire simulation process from continuous deformation to fracture motion.SDEM was jointly designed and developed by the computational element mechanics team of the Dalian University of Technology and the team at Dizaosoft,Dalian Co.,Ltd.It can effectively simulate the failure process in brittle materials such as sea ice and rocks (Ji et al.,2012).DICE2D is an open source discrete element software program developed by Professor Zhao Gaofeng of Tianjin University(Zhao,2015).

A high-performance discrete element software,MatDEM,was developed from scratch (Liu et al.,2013).By using an innovative matrix computing of the DEM,MatDEM realizes the high-performance discrete element numerical simulation of millions of elements in one computer.The basic principles of the DEM,the implementation and development of the MatDEM software,and its applications are introduced in this paper.

2.Basic principles of the DEM

2.1.Contact model of elements

As shown in Fig.1a,the DEM constructs a geotechnical model by depositing and cementing a series of elements with specific mechanical properties.The most basic linear elastic model assumes that the elements interact with each other via spring forces.The normal force (Fn) and normal deformation(Xn)between each two elements can be simulated by a normal spring(Place and Mora,1999):

whereKnis the normal stiffness of the spring,Xnis the normal relative displacement(Fig.1b),andXbis the fracture displacement.Initially,elements are interconnected with their adjacent elements and subjected to tensile or compressive spring forces(Equa.1a).WhenXnbetween the two elements exceeds the fracture displacement (Xb),the spring breaks and the inter-element tensile force no longer exist between them(Equa.1c);however,the compressive force may act between them when they return to a compressive status(Equa.1b).

Shear force (Fs) and shear deformation (Xs) between elements are simulated by tangential springs(Mora and Place,1993):

whereKsis the shear stiffness andXsis the relative shear displacement.

The spring also has a failure criterion in the tangential direction,which is based on the Moho–Coulomb criterion:

whereFsmaxis the maximum shear force,Fs0is the inter-element initial shear resistance,and μpis the inter-element coefficient of friction.

2.2.Damping force of elements

In order to simulate the process of energy dissipation in the real world and avoid the accumulation of stress wave energy in the system,the damping must be used in the DEM.Damping can be defined in several ways,among which the global damping(Fv)is the simplest one as follows(Liu et al.,2013):

wherex'is the element velocity at the current moment,and η is the viscosity.Therefore,the damping force on the element is equal to the damping coefficient times the velocity,and its direction is opposite to the velocity direction.This definition of damping force is adopted by the MatDEM.

2.3.Calculation of energy

The heat generated in the model includes viscous force,breaking energy,and sliding frictional heat,which correspond to damping,fracture,and friction,respectively.In a discrete element numerical simulation,the dynamic evolution of the discrete element system is simulated by time step iterations,and various kinds of heat are calculated and accumulated in each time step.

(1) Viscous force

When elastic waves propagate among elements,mechanical energy will gradually dissipate into heat due to friction,scattering,and other factors.In the DEM,damping is used to weaken the elastic wave energy of the model and dissipate the kinetic energy in the discrete element system.The damping force (Fd) acting on the elements is given by Equa.5:

Fig.1. Schematic diagram of the linear elastic model.

where η is the damping coefficient and ν is the element velocity.Because the time step of the simulation is relatively small,the element velocity in a step is assumed to be constant.Viscous heat can be calculated by Equa.6:

where dxis the displacement of the elements in the current time step.

(2) Breaking energy

Thus,the breaking heat generated is equal to the reduction of elastic potential energy.If the inter-element normal force is tensile,the elastic potential energy of both normal and shear springs will reduce to zero when the intact bond breaks,and the reducedEe,tis the sum of the elastic potential energy of both normal-and shear springs:

If the inter-element normal force is compressive,the normal spring force and the corresponding elastic potential energy will not be changed when the intact bond breaks.According to equation(3),the shear force reduced fromFsmaxtoFs’max,and the reducedEe,ccan be expressed as:

whereKsis the inter-element shear stiffness.

(3) Sliding Friction

Two elements begin sliding relative to one another when the bond is broken,and the magnitude of external force exceeds the maximum shear force of the broken bond.Frictional heat(Qf)generated during the sliding process is defined as the product of the average sliding friction and the effective sliding distance:

whereFs1andFs2are the sliding friction forces at respectively the beginning and the end of the current step(i.e.,the maximum shear forces of the broken bond),and dSis the effective sliding distance.Details about this equation can be found in) Liu et al.(2017).

3.The implementation of the MatDEM

3.1.Time stepping algorithm

Based on obtaining the force of each element,the displacement of the element is calculated using a time stepping algorithm.In a short step time(dT),the force,acceleration,velocity,and displacement of elements are calculated.Once the current time step calculation is completed,a time step is advanced to implement a new iteration of the DEM.The detailed calculation process is as follows:firstly,we calculate the current acceleration of elements via dividing the known resultant force by the mass based on the traditional Newtonian mechanics;then,in the current step time dT,we can obtain the initial velocity of the next time step by adding the current velocity plus the velocity increment(the acceleration multiplied by dT),and then calculate the corresponding element displacement;finally,we run a new iteration to realize the dynamic simulation of the DEM.

3.2.Differential calculation of the element motion

In the MatDEM,the linear elastic contact model is adopted by default,and global damping is adopted to consume energy.Essentially,the movement of elements in MatDEM can be summarized as the damping harmonic vibration of a spring oscillator.The equation for calculating the period(Tη)of the damping harmonic vibration is:

This formula contains two important clues:(1) the vibration period increases with an increasing element mass(m),a decreasing stiffness(k)or an increasing damping coefficient (η);(2) there is a critical damp-which produces an infinite period.When the damping coefficient is 0,we will obtain the equation of harmonic vibration

Fig.2 shows the displacement-time curve for one period of the harmonic vibration of an element.The numerical simulation of DEM is based on time stepping iteration,which assumes that the acceleration of the element is constant within a small time step.Therefore,it is necessary to select a time step that is short enough to accurately simulate a harmonic vibration process.To properly fit a sinusoidal displacement curve such as Fig.2,the time step must be lower than 10% of the period.A series of numerical simulation results indicate that when the time step is 2%of the period,the regular calculation speed and accuracy can be satisfied,which is used in the MatDEM by default.

3.3.Element damping and optimal damping coefficient

A sample codeTwoBallsof the MatDEM can be used to validate the discrete element simulation.This code simulates the process of a damping harmonic vibration.Fig.3 shows the simulation results for two identical elements vibrating on fixed boundaries.Fig.3b displays the energy curves for simulation results with the critical damping,which allows each element to return to equilibrium as quickly as possible without periodic vibration.However,if 1/10 of the critical damping is adopted (Fig.3c),the kinetic energy will decrease gradually during the fluctuation.

Fig.2. Element motion during a vibration period.

Fig.3. The damping harmonic vibration process of an element.(a)The model;(b)Energy curves when using the critical damping;(c)Energy curves when using 10%critical damping.

In the discrete element simulation of quasi-static processes (such compressive strength tests),a greater damping coefficient is usually adopted to quickly consume the kinetic energy of the discrete element model,thereby reducing the computational time of the numerical simulation.In numerical simulations of tests,the loading speed can be simulated by adjusting the damping coefficient and the iteration step number.If a displacement is applied to an element with a critical damping coefficient of 1/10 and calculated iteratively for 100 times with a new displacement load added immediately before the stress wave dissipates,the compression process is considered a fast loading (equivalent to the“fast”in Fig.3c).If a displacement is applied to an element and then run 1000 times of iterations,the stress wave and kinetic energy will be dissipated and a slow compression or quasi-static effect (equivalent to the“slow”of Fig.3c) can be achieved.Moreover,if the critical damping is adopted,quasi-static simulation can be achieved after 100 iterations(Fig.3b).

Based on the critical damping formula,an optimal damping for multielement systems can be obtained using the following semi-empirical formula(Liu et al.,2017a):

whereVis the model volume;dis the element diameter;mis the element mass;andkis the element stiffness.In quasi-static issues(such as a test of uniaxial compressive strength),the optimal damping coefficient is adopted to rapidly converge the system energy,which equals the critical damping of an element divided by the number of elements in one direction of the model.

3.4.High-performance matrix computing of the discrete element method

The DEM decomposes rock and soil mass into a series of elements spatially,and performs iterative calculation in time;such method is computationally intensive.The number of calculation elements for a three-dimensional numerical simulation is usually limited to tens of thousands,resulting in the sample size being too small to meet the actual engineering requirements.

Recent studies have shown that high-performance computing technology based on a graphics processing unit(GPU)can effectively improve the speed of numerical calculation.It can improve computational efficiency by 10–50 times and has been widely used in molecular dynamics,petroleum exploration,computational mechanics,materials science,and artificial intelligence.The GPU is commonly referred to as a graphics card,such as NVIDIA's discrete graphics card that contains the CUDA computing cores needed for the GPU computing.Ordinary graphics cards usually have hundreds of CUDA computing cores.In particular,professional Tesla GPUs,such as P100 and V100,are equipped with thousands of calculation cores that vastly improve the computational efficiency.In view of the immense computational complexity of the DEM,the authors independently developed the high-performance discrete element software MatDEM.The discrete element simulation of millions of elements was realized by adopting the innovative matrix discrete element computing method and a three-dimensional contact detection algorithm.This enhances the calculation efficiency and the number of elements reaches 30 times of that of the commercial software,and as a result it can be used in the engineering applications.

Fig.4. Calculation speed of CPU and GPU with different numbers of elements.

Fig.4a shows the curves of tested speeds of CPU and GPU on an ordinary laptop (ThinkPad T470).This shows that when the number of elements is small,the GPU calculation speed is lower than that of the CPU because most of the calculation time is spent on communication between them.As the number of elements increases,the GPU computing speed increases rapidly.When the number of elements exceeds 10 000,the calculation speed of the GPU is significantly higher than that of the CPU.When the number is 100 000,the calculation speed of the GPU is about 1 million element movements per second,five times that of the CPU.Fig.4b shows the speed test results for the MatDEM on a professional GPU server.Similarly,when the number of elements exceeds 10 000,the GPU speed is significantly higher than that of the CPU,and when the number of elements reaches 1 million,the GPU calculation speed is about 50 times that of the CPU.As shown in Fig.4b,when the number of elements exceeds 1 000,the CPU calculation speed is basically unchanged while the GPU speed increases rapidly as the number of elements increases.Therefore,the GPU matrix computation is especially suitable for large-scale discrete element simulation.

The calculation speed of the MatDEM increases rapidly when using the calculation method for the GPU matrix in line with the increase in the number of elements.Numerical simulations of the DEM are based on repeated iterative calculation,and the number of iterations per second determines the time of calculation.A professional GPU server (Tesla P100 GPU,16 GB memory)supports 1.5 million 3D elements at most,and there will be a high number of iterations per second(more than 52 times)within 200 000 elements.For a high-performance desktop computer(GTX 1080Ti),there will be a high number of iterations per second(29 or more) within 100 000 elements.For regular laptops (GT 940M),there will be a high number of iterations per second (more than 25 times)within 30 000 elements.Table 1 gives the optimal interval of element numbers for these video cards,which suggests that the number of computing elements in the test can be calculated more quickly.

Table 1 also shows the maximum number of simulation elements in each GPU card(or video card)and the recommended number of elements of high-precision simulation.The MatDEM needs approximately 1 GB memory per 100 000 elements of simulation.The latest Tesla A100(48 GB memory)currently supports approximately 4.5 million elements.

Table 1 Number of elements simulated by the MatDEM when using different GPU cards.

4.The software MatDEM and its applications

4.1.Brief introduction of the MatDEM

The MatDEM has realized the automatic modeling of discrete element materials,multi-field and fluid-solid coupling numerical simulations,the energy conservation calculation of discrete element system,etc.It integrates pre-processing,high-performance solver,post-processing,and secondary development to facilitate the large-scale discrete element simulation of various geological and engineering problems.

Fig.5 is the startup form of the MatDEM 2.0.The software supports 6 different languages.The major functions and characteristics of the Mat-DEM are listed as follows.(1)Preprocessor:Support reading and writing of multiple types of files,digital elevation and image modeling,provide rich geometric modeling functions,etc.(2)Solver:Support both GPU and GPU computing;support 3 million elements in one GPU.(3) Post processor:provide dozens of drawings,the results can be saved as images and vector files,GIF animations and videos.(4) Extension module:powerful scalability,self-service addition of new attributes,pre-andpost-processing functions,functions and window applications,support for custom discrete element contact models.(5)Secondary development:Provide dozens of case codes and App source code,including test simulations,geological disasters,geotechnical engineering,dynamic actions,planetary motion,multi-field coupling,fluid-solid coupling and other categories.(6) Real-time control:Users can stop the calculation at any time,view all the calculation parameters,and then continue the calculation process.The software provides over 200 APIs (Application Programming Interface) and over 40 sample source codes.The number of lines of these code is usually between 100 and 300 lines.By using the sample codes,one may build numerical models and simulate various processes.

4.2.Applications of the MatDEM

MatDEM has been successfully applied to a range of engineering and geological studies,including the formation of compaction bands (Liu et al.,2015),the analysis of foundation pit stability (Suo et al.,2017),tunnel issues (e.g.;Xue et al.,2020),confined compression test simulations of sand (e.g.Qin et al.,2018;Zhu et al.,2020),landslide simulations(e.g.Scaringi et al.,2018;Liang et al.,2019;Zhu et al.,2020),etc.Fig.6 shows the simulation results of the four typical cases,including seismic wave propagation across a V-shaped valley,hydraulic fracturing,rock breaking by TBM cutter and landslides.

A V-shaped valley can produce ground motion amplification effects,and aggravate the destructive disaster effects of earthquakes.The case shown in Fig.6a simulates the propagation of seismic wave in a valley model that contains 107 thousand elements.A compressive wave was generated on the left boundary of the model,and propagate across the valley.The simulation results show that the back wave slope produces an obvious amplification effect of ground motion acceleration,which is extremely intense near the bottom side of the valley.

Fig.6b shows the numerical results of hydraulic fracturing.The material used in the model is sandstone.Four boundaries of the model are fixed,and a constant pressure of 100 MPa is applied at the center of the model.Three micro cracks are generated after 0.02 s,and the angle between them is approximately 120°(Liu et al.,2020).The color in the figure indicates the pore pressure.

The MatDEM also has been applied in the numerical simulation of rock breaking by a TBM cutter(Liang et al.,2019).In Fig.6c,a cutter and a soft-strong inter-bedded rock model with 502 thousand elements was built.In the simulation,the cutter rolls on the rock block surface from left side to the right side.The forces on the cutter,element connections,energy conversion and heat generation are monitored,providing important clues for deciphering the dynamic process.

As shown in Fig.6d,the software is used to construct a threedimensional discrete element model based on digital elevation data,and simulate and analyze the movement and accumulation process of a landslide.The total number of elements in the landslide model is 826 000.The numerical simulation includes 432 500 steps,which simulate the real time for 155.7 s.By using a GPU workstation(Tesla P100 GPU)for the calculation,this simulation takes about 18 hr.

5.Development and prospect of the MatDEM

The author began the investigation of the DEM algorithm during postgraduate study at Nanjing University in 2007,and developed the original theoretical test software using the C#language.In 2011,based on the idea of the matrix computing,a solver of 2D discrete element method was developed.After that,the development of 3D discrete element software began in 2012.In 2013,the high-performance matrix computing of the DEM was realized,which greatly improved the computing efficiency by dozens of times.Furthermore,the number of computing elements reached the level of million.

Fig.5. Startup form of the MatDEM 2.0.

Fig.6. (a) Seismic wave propagation across a V-shaped valley;(b) hydraulic fracturing;(c) rock breaking by a TBM cutter;(d) landslides.

From 2014 to 2016,the MatDEM has realized a large-scale discrete element heat calculation and energy conservation calculation,and established a complete post-processing module as commercial software.The multi-field coupling simulation has gradually being improved and has been applied to the simulations of geological and rock mechanics issues.Since 2017,with the support of the Chinese Society of Rock Mechanics and Engineering,the software is used in engineering applications.On May 4th-6th,2018,during the forum jointly hosted by Chinese Society of Rock Mechanics and Engineering and Nanjing University,the MatDEM version 1.0 was officially released.In 2019,the software received“The Digital Simulation Independent Software Innovation Award of China”.Until November 2020,there are over 10 thousand downloads of the MatDEM from over 30 countries.Over one thousand users are using the MatDEM in their studies and projects.

Recently,a pore density flow method is proposed and is included in the MatDEM,which may simulate the fluid flow at a pore scale,and the interaction between fluid and solid phases (Liu et al.,2020).It may be used to simulate many complicated physical processes and geological phenomena,such as magma intrusion,hydraulic fracturing,and engineering grouting,as well as promoting the development of research in numerous fields.

The software,source codes and tutorials (in both Chinese and English) are provided online http://matdem.com.In addition,there are textbooks in Chinese (Liu et al.,2017b) and English (published in 2021 by Springer).The MatDEM team will continue to optimize the MatDEM computing core and provide new standards and open API.The team are also pleasured to work with experts in various fields,and develop international high-performance discrete element software for both the academic community and industry.

Acknowledgments

Financial supports from the Natural Science Foundation of China(41761134089,41977218),Six Talent Peaks Project of Jiangsu Province(RJFW-003),and the Fundamental Research Funds for the Central Universities (14380103) are gratefully acknowledged.We are grateful to anonymous reviewers for many valuable comments that notably improved the manuscript.