A robust potential-based contact force solution approach for discontinuous deformation analysis of irregular convex polygonal block/particle systems

Contact interaction of two bodies can be modeled using the penalty function approach while its accuracy and robustness are directly associated with the geometry of contact bodies. Particularly, in the research fields of rock mechanics, we need to treat polygonal shapes such as mineral grains/particles at a mesoscale and rock blocks at a macroscale. The irregular shapes (e.g., polygons with small angles or small edges) pose challenges to traditional contact solution approach in terms of algorithmic robustness and complexity. This paper proposed a robust potential-based penalty function approach to solve contact of polygonal particles/block. An improved potential function is proposed considering irregular polygonal shapes. A contact detection procedure based on the entrance block concept is presented, followed by a numerical integral algorithm to compute the contact force. The proposed contact detection approach is implemented into discontinuous deformation analysis with an explicit formulation. The accuracy and robustness of the proposed contact detection approach are verified by benchmarking examples. The potential of the proposed approach in analysis of kinetic behavior of complex polygonal block systems is shown by two application examples. It can be applied in any discontinuous computation models using stepwise contact force-based solution procedures.

Block rotation angle about axed oz e x ; e y ; c xy Constant normal and shear strain of a block T i x p ; y p À Á 2 Â 6 transform matrix of point p in block i t n Beginning of step n t nþ1 Ending of step n t nþ1=2 Half step instant of step n t nÀ1=2 Half step instant of step n À

Introduction
Contact occurs almost everywhere in our daily life [29,42], and it is involved in scientific and engineering problems with a broad scale range (e.g., interaction of particle in air with a scale of lm, collision of planet and comet in space with a scale of km). Particularly, contact of mesoscale mineral particles and macroscale rock blocks is of great important in the research fields of rock mechanics, as the strong relevance of geometrical and physical properties of the discrete bodies and their interfaces to the deformation and strength of rock material and the safety of engineering activities. In the research of the relationship of mesoscale grain properties and macroscale rock deformation and strength properties, the Voronoi-based element shape is commonly used [19,38,46,53]. In the analysis of rock block systems, an assembly of polygonal block is often formed by the intersection of joints/boundaries [9,11,20,34,45,54,58]. Compared to circular/spherical shape, the polygonal/polyhedral shape brings some difficulties for mechanical analysis of block/particle assembles, such as treatment of non-smooth contact normal change at vertices [1,5,40,55,61,62], treatment of quasi-parallel edges [63], and the contact type detection of irregular polygons with small edges or small angles.
Solution of contact should consider both the geometrical and physical properties of the contact interfaces. The contact state and contact force are determined by the geometrical relationship of two bodies and the interface constitutive model, which lead to strong nonlinear properties. Many numerical methods have been proposed to solve the mechanical and coupled hydro-thermal-mechanical response of a model with discontinuities by indirect [37,66,67] or direct consideration of contact [21,22,43,50]. Among them, the ''discrete element'' approaches allow both finite displacement, rotation, and detachment of discrete bodies and recognize new contacts automatically [3]. The commonly used approaches belonging to the discontinuous type include rigid block spring method (RBSM) [52], discrete element method (DEM) [2], discontinuous deformation analysis (DDA) [6,24,33,59], the combined finite-discrete element method (FDEM) [26], and the numerical manifold method (NMM) [13,18,23]. The solution strategies of these discontinuous models can be mainly classified into two types: the explicit contact force-based solution approach and the implicit displacement-based solution approach. The explicit approaches detect contacts and compute contact forces for the integration of block/particle displacement in each step. By contrast, the implicit approaches detect contacts and solve a global equation of motion that consider the inequality contact constraints. The explicit approaches are restricted by a critical time step, while the implicit approaches require the convergence of contact states. In this paper, the focus is given on the accuracy and efficiency of the explicit approach in solution of polygonal block/particle systems with irregular shapes.
The major issues in explicit contact analysis approach are the detection of contact and computation of contact force. The contact detection is usually proceeded using a two-phase strategy that includes a rough search phase to locate neighbor block pairs and a delicate search phase to find contact types or obtain overlapping volumes [60]. The algorithms to detect neighbor block pairs are well established, such as the no-binary search [27], double-ended spatial sorting [31], and CGrid algorithms [41]. In this paper, the focus is given on the delicate search algorithms for polygons, which can be classified into two types: (1) type-based search and (2) overlap-based search. As stated in [1,5,55,63], there are some embedded difficulties in treating singularity at vertices and irregular shapes. By contrast, the overlap-based detection algorithm aims to establish a structure of overlapping area (in two-dimensional problems) or overlapping volume (in three-dimensional problems), and contact force can be computed using the penalty function approach on the basis of potential [8,12,28] of the contact bodies. In this paper, the overlapbased contact solution approach is used in view of its advantages in treating polygonal shape-induced difficulties.
The potential-based penalty function approach [28] provides a robust framework for contact force computation of two contacting bodies. Some recent improvements include the redefined potential function for triangles with clear physical meaning [49], the improved distance potential for convex polygons [56], concave polygons [7] and convex polyhedra [57], the potential traction based on triangle meshes [47] and tetrahedral element [48], robust potential function for irregular polyhedra [65], and the smooth contact algorithm [17]. In this paper, the original contact potential-based penalty function method [28] is improved to be suitable for polygons with small edges, small angles, or small faces. The improvement originates from the basic idea ''a penetrated point in a polygon will be pushed out along the shortest path.'' On this basis, a new contact potential function is defined and a novel contact force computation approach is proposed, including the geometrical associated algorithms of overlap examination and intersection polygon construction, and the physical associated approach to compute the contact force with a numerical integral scheme. This paper begins with the formulation of the explicit DDA approach based on explicit representation of contact force. Then, the improvements on the potential function methods and the computation of contact force (normal forces and tangential forces) are introduced, with particular focus on some key algorithms in contact detection and contact force computation process. The accuracy and robustness of the improved approach are verified by benchmarking contact scenarios, and its potential in kinetic analysis of complex block systems is further investigated by the application examples. A summary of the proposed approach is given in the last section.
2 Explicit solution approach based on potential contact force The contact force-based solution approach is applied in the original framework of two-dimensional discontinuous deformation analysis. The displacement function, equation of motion, and stepwise solution procedure with the explicit formulation are briefly introduced.

Displacement function
The linear displacement assumption in original DDA program is still applied in this paper. The displacement u p x p ; y p À Á of point p in block i is computed by where T i x p ; y p À Á is a transform matrix similar to the shape function matrix with dimensions of 2 Â 6 and d i is the degree of freedom in terms of block centroid displacement and block constant strain. d i is represented by where u o ; v o ; r o are block centroid translation along axes ox and oy; r o is rotation angle about axes oz; and e x ; e y and c xy are the constant normal and shear strain. The transform matrix T i x p ; y p À Á is computed by where x t ¼ x p À x o , y t ¼ y p À y o , and x o ; y o ð Þ denotes the centroid coordinate of block i.

Equation of motion
Kinematics and dynamics of a block/particle system are characterized by the contact interaction mechanism of discrete bodies. The kinetic properties of a deformable block/particle can be described by the equation of motion, which can be established by minimizing the potential energy [4,33], where M i , C i , K i , and F i represent the mass matrix, damping matrix, stiffness matrix, and force vectors of block i; the dimension of the matrix and vectors are 6 Â 6 and 6 Â 1 respectively; and _ d i , € d i are the first-order and second-order derivative of the unknown vector d i .
M i , C i , and K i can be obtained by the following integral format, where T i is the transform matrix as in Eq. (3); E is the elastic stiffness matrix; and a and b are two coefficients associated with the damping proportional to the mass and stiff terms. F i includes all external forces (e.g., point force, surface stress), contact force vectors, and initial stress force vectors. Details of the above matrix can refer to [33].

An explicit approach
The global equation of motion of a block system can be established by minimizing the global potential energy function that includes the contact potential energy with the penalty function format [64]. The original DDA uses an implicit time integration approach to solve the global equation of motion, which needs to solve the global matrix each step [15]. In addition, some explicit time integration approaches that directly use contact forces to solve Eq. (4) for each block separately are also proposed [7,25,30,47,51,64]. For the computation of incremental displacement from time t n to t nþ1 , the velocity Verlet approach [51] solves Eq. (4) at t n and use the obtained acceleration term to compute the incremental displacement; the predictor-corrector approach and its modified version [64] solve Eq. (4) at t nþ1 with a predicted contact force and then correct the incremental displacement; and the direct approach solves Eq. (4) at t nþ1 with the contact force obtained at t n using the either the Newmark integration approach [30,47] or generalized-a approach [7]. A comprehensive comparison of the velocity Verlet approach, constant acceleration integration approach, and the predictor-corrector solution approach can be found in [16,64]. In this paper, the velocity Verlet approach [32,51] is used for time integration of the equation of motion. For the analysis of the physical process from time t n to where K i d n i in Eq. (4) is added to F n i with the vector format of step initial stress. In this paper, C i _ d n i is not considered and details of the coefficient choice can refer to [14]. The acceleration term, € d n i , can be obtained by solving Eq. (6); the velocity at t nþ1=2 ¼ t n þ h=2 ð Þcan be can be computed by The incremental displacement in this step can be obtained by Note that the model configuration is updated each step and d n i in Eq. (8) becomes zero automatically. Given the Initial condition and boundary conditions, the deformation and displacement of a discrete body system can be obtained stepwise using Eqs. (6)(7)(8). In this paper, the focus is given on the formulation of contact force and its contribution to force term F n i in Eq. (6).

The improved potential-based contact approach
The concept of ''contact potential'' is initially proposed by [28] in their work on penalty function method for combined finite-discrete element method, and it is used for contact force computation among triangles and tetrahedrons. The definition of potential has been improved for triangles [49], convex polygons [56], convex polyhedra [57], concave polygons [7], and irregular polyhedra [65]. The finite element topology was also applied to construct a smooth potential field in FDEM simulations [17]. However, it still lacks a unified and robust definition for polygonal shapes when irregular shapes (such as small angles and small edges) are encountered. To overcome these shortcomings, an improved potential definition is proposed and the work is detailed in this section. The contact solution procedure in explicit solution approach can be conducted with two steps: (1) a contact detection step to judge the body overlap and to construct the intersection shape, and (2) a force computation step to obtain the contact force based on the intersection shape and the body potential. The improvement toward these two steps is detailed in Sects. 3.1 and 3.2, respectively. Moreover, a solution procedure for contacts of polygons and boundaries is discussed in Sect. 3.3.

Detection of contact
The first step of a contact solution procedure is the judgment of contact states and extraction of the geometrical information. For the proposed model, it consists of a threesub-step procedure: (1) neighbor search; (2) overlap examination; and (3) intersection polyhedron construction.

Neighbor search
This is a common phase for all kinds of approaches that need to solve contacts of discrete body assembly. In this phase, a neighbor body list will be established for each discrete body. The efficiency of this detection phase becomes significantly important with increasing body numbers. Many algorithms have been published to enhance the computational efficiency in neighbor search process, such as the grid-based searching algorithm [2], the no-binary search algorithms [27] for discrete bodies of similar sizes, the double-ended spatial sorting algorithm [31], the CGrid algorithm [41] for discrete bodies of various sizes, and the multi-shell cover algorithm for polyhedra [39,44]. The neighbor search algorithm should be chosen according to the distribution feature of the discrete body system.

Overlap examination
Although some polygons are marked as neighbor pairs after detection of neighbor block pair, they may not contact at the current configuration. Overlap examination is a detection step to exclude the neighbor polygon pairs that are near but no-in-contact. The next step to compute the intersection polygon will only be triggered if overlap occurs for a neighbor polygon pair.
Overlap checking of two close bodies is a typical topic in computational graphics and discontinuous computation methods. It is essentially, from a mathematical view, relationship judgment of two sets. This judgment can be simplified with the concept of entrance blocks [35]. The overlap judgment of two polygons can be simplified into the examination of a chosen reference point and their entrance block. For two polygons termed as polygon A and polygon B, choosing the centroid of polygon A as the reference point a o , the entrance block E a; b ð Þ is defined by It has been proofed that the entrance block of two convex polygons is still convex polygons [35]. Overlap status of two convex polygons A and B can be judged by a point-in-polygon examination using reference point a 0 and entrance block E A; B ð Þ. The overlap examination is proceeded in two steps. The first step is to obtain the boundaries of the entrance block that formed by an entrance block computation of the vertices and edges from A and B, respectively. In this step, some vertex-edge pairs that do not form valid entrance candidate [61] will be omitted. Algorithm detail of this step is shown in  Table 2. An example to check overlap of two convex polygons is shown in Fig. 1. In this example, block B is fixed, and position of block A 2 is obtained by translating A 1 .
Þ when the centroids c 1 and c 2 of polygon A 1 and A 1 are used as the reference point. As illustrated in Fig. 1,

Intersection computation
The intersection polygon (IP) will be constructed if overlap of two polygons occurs. As it is prior known that the intersection polygon of two convex polygons is still Table 1 The proposed algorithm to form entrance block boundaries of two convex polygons Table 2 The proposed algorithm to judge position of reference point a 0 and convex entrance block E A; B ð Þ Fig. 1 Illustrations of overlap judgment based on entrance block concept convex, a simple polygon cutting algorithm can be used to construct the intersection polygon. The two blocks from the neighbor block pair are termed as contact polygon (CP) and target polygon (TP), respectively. The lines extended from each edges of CP will be cycled to cut the TP, and in each cycle the newly generated polygon that locates at the half space below the cutting line will be stored for further operation. This cutting process is illustrated in Fig. 2. After the cycling of each CP edges, the remaining polygon is IP. Details of the algorithm are shown in Table 3.

Potential contact force
For a neighbor polygon pair, choosing anyone as the contact polygon (CP), the remaining one is the target polygon (TP). If overlap of CP and TP occurs, their intersection forms the intersection polygon (IP). As shown in Fig. 3, the global normal contact force f ng CP acting on CP can be computed by integrating df n CP on IP, where u CP and u TP are the potential functions for the CP and TP, respectively; S IP denotes the area of IP. The global contact force f ng TP on TP has the same magnitude and opposite direction compared to f ng CP . As u CP and u TP are continuous functions, the contact force f ng CP can be equivalently computed by the integral over the IP boundary L IP [28], where n IP denotes the unit normal vector of a boundary integral point that point outside of IP.

New potential function definition
In analysis of a rock block system or a mineral particle assembly, there might be irregular polygons with edges or angles smaller than a characteristic value. These irregular shapes can cause inaccuracy, inefficiency, or errors in some existing contact solution approaches. Specially, for the potential-based contact solution approaches, there still lack a robust potential function definition for arbitrarily shaped polygons with above-mentioned irregular shapes. A new definition of potential function is proposed in this section, and the results are compared with previous approaches [56].
In this paper, the potential function for a polygon is redefined based on the concept ''a point penetrated in a polygon will always be pushed out along the shortest path.'' For a point belonging to a polygon, once the shortest distance d exit to exit a polygon is determined, the potential at that position can be obtained. The potential of a point p in a polygon X b is defined as where k n is the normal penalty parameter; H is a unit length value that can be set to average length of all polygon edges of the assembly; d exit is the shortest distance for point p to exit polygon X b . To keep energy conservation in contact, the potential values of any points on polygon boundary should be equal and the definition satisfies this requirement. Equation (12) gives the potential function definition of both convex and concave polygons. As the potential function values in IP will be used for computation of contact force, the difficulty lies in how to obtain the distribution of the potential function in an arbitrarily shaped polygon. Munjiza and Andrew [28], Yan and Zheng [49], and Xu et al. [47] show the computation of potential function values in a triangle. Zhao et al. [56] give the computation methods for potential function of convex polygon based on the subdivision strategy. Fan et al. [7] give the computation of potential function for arbitrarily shaped polygon based on the classification of contact type. However, these definition and computation approach may lead to robustness problem or inaccuracy when polygon with small angles or small edges is involved in the block assembly.
Using this definition for convex polygons, the shortest distance d st can be obtained by where e i represent boundary edge i of the polygon X b ; p is the point whose potential value is to be computed; and d e i p The potential of a point in a convex polygon can be computed conveniently using the proposed definition in Eqs. (12)- (14), even for irregular shapes with small edges or small angles. Additionally, with the same H value in the potential computation of each polygon, it ensures the contact potential of the points with the same penetrations keeps equal. Figures 4 and 5 show some examples of the potential distribution in polygons. In these examples, k n ¼ 1 Pa and H ¼ 1 m. Figure 4 shows the comparison of the potential distribution in a polygon with small edges with the proposed function (12) and the previous definition [56]. For a polygon with small edges, there would be a large inner part without a valid potential definition when the previous definition is used. This problem may become significant when relatively large penetration is allowed in analysis of a block system with irregular polygonal shapes. By contrast, the proposed definition ensures the potential value can be obtained for each point inside a polygon.     Fig. 5, the potential distribution in a concave polygon using the definition of (12) is shown. This distribution is obtained by decomposing the concave polygon into six convex ones by demarcation lines. However, this paper focuses on convex shapes and decomposition of an arbitrarily shaped concave polygon is beyond our scope.

Numerical integral segments
With the definition in Eq. (12), the distribution of potential value on boundary edges of IP is piecewise linear. The key issue for the integral of Eq. (11) is to find all demarcation points which divide the potential field along an edge into several linear distribution pieces. If the all demarcation points on an edge is obtained, the integral of potential on an edge can be computed accurately piece by piece. Then, the total contact force can be obtained by adding the accumulated force vectors of all IP edges.
The effects of small edges are considered in this paper. If the length of a CP or TP edge is smaller than a specified tolerance, a virtual polygon node will be generated by intersecting the two bisection projection lines of the two angles neighboring the small edge. Both the real nodes and virtual nodes of CP or TP are used in locating the demarcation points.
The potential demarcation points on a boundary edge of IP are obtained as follows:  (4); if the centroid point is not on boundary of TP, the intersection points of the edge and the rays of angle bisector starting from TP nodes and TP virtual nodes will be added. (4) All points obtained from steps (1) to (3) will be sorted from the beginning to the ending of the IP edge. Figure 6 illustrates the demarcation points computation process of an IP edge corresponding to a polygon (either CP or TP). If the polygon is CP, the demarcation points include beginning, central, and ending points of the edge, and all intersection points of the edge and rays of angle bisector starting from nodes or virtual nodes of CP. After locating all demarcation points, the edge will be naturally divided into several integral segments.

Formulation of normal force
If boundary edge e i of IP are discretized into sn segments with index sk ¼ 1; 2; . . .; sn, the integral of normal contact force on edge e i can be computed by Assume points p a and p b are the two ending points of segment sk. As the potential is linearly distributed along segment sk, the integral of potential force on sk can be computed by where f a ¼ u CP p a ð Þ À u TP p a ð Þ; f b ¼ u CP p b ð Þ À u TP p b ð Þ; L sk is the length of segment sk. Meanwhile, the coordinate vector p k ¼ x k ; y k ð Þ of the loading point p k on segment sk can be obtained by where p a and p b are the coordinate vector of point p a and p b , respectively. The distributed normal contact force can be equivalently represented by a total normal contact force f ng CP , where ei indexing the contact force f n ei of edge e i and en is the total number of IP boundary edges; sk indexing the contact force f n sk of segment sk and sn is the total number of segments on edge e i . The unit contact normal vector n g of f ng CP can be obtained by If p m ¼ x m ; y m ð Þ is a point on the acting line l ng along which the moment of total normal contact force is zero, it satisfies the moment balance equation along axes oz, where f nx k , f ny k represent components along ox and oy direction of the normal force on segment sk; snn is the total segment number of the IP boundary. x k ; y k ð Þ is the coordinate of contact force f n sk . The loading point p g ¼ x g ; y g À Á of the global normal contact force is obtained as follows: (1) an auxiliary line l tg that passes the IP centroid and perpendicular to acting line l ng is created; (2) the intersection point of the auxiliary line l tg and the acting line l ng is chosen as the contact point of global normal contact force and global tangential contact force. The normal contact force works along l ng , while the tangential contact force works along line l tg .
Recall the formulation of equation of motion using the potential energy format, and based on the block frame at the beginning of the step, the potential energy of the contact force for CP can be represented by For CP, the contribution of the contact force to the force terms F n i in Eq. (6) is Similarly, for TP, the contribution of contact force f ng TP to F n i is

Formulation of tangential force
The global tangential contact force is represented in a global format, and Coulomb friction law is used to differentiate the stick and sliding states of two contacting polygon. For a block pair CP and TP, the accumulated relative displacement d will be recorded. Stick status is assumed if they come into contact initially in a step and d at that instant is set to zero. The normal and tangential components of d are named d nor and d tan , and the relative sliding direction is s. At instant t n , for a contact pair with stick status in previous step (from t nÀ1 to t n ), d n nor , d n tan , and s n are determined by d n nor ¼ d n n gn ð Þn gn ð24aÞ The tangential contact force f tg CP for CP is obtained by when kk tan d n tan k\kf ng CP k tan / Àkf ng CP k tan /s n when kk tan d n tan k ! kf ng where / represents the frictional angle; k tan is the tangential penalty value. At instant t n , for a contact pair with sliding status in previous step (from t nÀ1 to t n ), Eq. (24) is still applied for computation of d n nor , d n tan , and s n . Assume the frictional force for CP at t nÀ1 is f tg nÀ1 , the tangential contact force f tg CP for CP at t n is obtained by Based on the frame at t n , the potential energy contributed by tangential contact force for CP is For CP, The contribution of tangent contact force f tg CP to the global force F n i is For TP, the tangential force f tg TP ¼ Àf tg CP . Similarly, the contribution of tangential contact force f tg TP to global force F n i is If the contact status at instant t is sliding, d t will also be set to zero.

Contact information update
For process from t n to t nþ1 , once the incremental displacement d i for each block is obtained, the incremental displacement Du g CP and Du g TP for contact point p g in CP and TP can be computed, and the incremental relative displacement Dd n ¼ Du g CP À Du g TP . The accumulation of relative movement can also be updated. For stick contact pairs at t n , the accumulated relative displacement at t nþ1 d nþ1 ¼ d n þ Dd n ; for sliding contact pairs at t n , the accumulated relative displacement at t nþ1 d nþ1 ¼ Dd n .

Polygon-boundary contact
Contact of polygon and model boundary needs to be considered in some simulations. Two kinds of boundary conditions are used in current models. The first kind of boundary is the usage of virtual polygons. These virtual polygons are involved in contact detection, while their displacement and deformation are not updated. The second kind is boundary lines and segments. Similarly, these lines and segments are regarded as parts of some virtual polygon boundaries and contact of polygon and lines are treated as contacts of polygons.

Verification of collision and frictional problems
In this section, the accuracy of the proposed contact solution approach for collision and frictional sliding problems is verified.

Collision of two identical polygons
In this example, collision process of two identical quadrilateral blocks was simulated using the proposed approach. Two square blocks with an edge length of 1 m are established, and their centroid is in a same line along ox axes, as shown in Fig. 7. The block on the left side is named block A, while the block on the right side is named block B. Centroid of block A has an initial velocity of 1 m/s along ox direction, while block B is static. The parameters for block physics and contact are shown in Table 4. The velocity of block centroid along ox direction is tracked, and the velocity change with time is shown in Fig. 8. It can be seen that the two blocks exchange their kinetic energy during the collision process and the results satisfy the conservation of momentum.

Collision of a polygon and a symmetric polygonal system
In this example, the collision process of a block and a symmetric block system was investigated with the proposed approach. Figure 9 shows details of the model geometry. The model geometry is symmetric regarding the oy axes passing the centroid of the top block. It can be predicted the block system will move in a symmetric way during the collision process. The parameters used in simulation are shown in Table 4. The model geometry and block velocity during the collision process are shown in Fig. 10. The movements of all block centroids were tracked during the simulation, and their traces are shown in Fig. 10. The symmetric response of the block system verifies the correctness of the proposed contact solution approach.

Sliding block on an inclined plane
In this example, the accuracy of the proposed algorithm for computation of block sliding displacement on an inclined plane is investigated. The model is shown in Fig. 11. Block B can slide along the inclined plane of block A if the friction angle / is smaller than the inclined angle h. If the initial velocity of B is zero, its sliding displacement l along the inclined plane can be computed analytically by This example is computed by the proposed approach with the parameters shown in Table 4. The values of friction angle / ¼ 0 ; 10 ; 20 are used to test the accuracy of the numerical approach. The time-displacement curves under various friction angle are shown in Fig. 12. Good agreements are obtained for numerical results and analytical results.

Robustness test examples
In this section, the robustness of the proposed contact solution approach in kinetic analysis of block system with various polyhedral shapes is tested and analyzed by several benchmarking examples. The physical parameters in simulation examples in Sect. 5 are shown in Table 5.

Collision of polygons with different sizes
In this example, the bouncing block model with various block sizes was tested to show the necessity of using unified H value in potential definition (12). Two cases are designed with different fixed block sizes, as shown in      8 Velocity changes of two blocks along ox direction during the simulated collision process Fig. 9 Initial geometry of the symmetric polygonal system significantly different block bouncing response for these two testing cases.

Polygon with small angles
In this example, collision of polygons with small angles was tested. Two identical triangles with a small angle of 5 are generated, as shown in Fig. 15. The two triangles move upwards and downwards, respectively, with an initial velocity of 1 m/s. The collision process is modeled by the proposed approach with the parameters shown in Table 5.
The two angle end points and polygon centroid are tracked to show the collision process, and the final configuration after 0.2 s is shown in Fig. 16. The result shows the two triangles move symmetrically regarding the collision centroid, verifying both the accuracy and robustness of the proposed approach in treating small angles. Fig. 13 The two bouncing block cases with different fixed block sizes

Polygon with small edges
In this example, collision of polygons with small edges was tested. Two polygons are generated and the bottom block has a small edge of 0.01 m compared to the typical block size scale around 1 m. For this case, when block penetrations around 0.01 m is allowed in the computation, part of the intersection polygon may locate to the inner parts without potential function definition when previous approach [56] is applied. This will cause the routine errors or inaccuracy of the computation. By contrast, these robustness issues are eliminated using the proposed approach with new potential function definition.
In this example, two model configurations are considered, as shown in Fig. 17. For the configuration in Fig. 17a, the vertex from the upper block will touch the small edge. For the configuration in Fig. 17b, the vertex of the upper block will touch the right edges neighboring the small edge. The physical parameters in simulations are shown in Table 5. The configurations after 0.2 s are shown in Fig. 18a, b for the two cases in Fig. 17a, b. The tracing of contacting vertex is also shown in Fig. 18. The results show the proposed approach can treat contacts of small edges with good robustness.

Application examples
The stability of rock structures on the earth surfaces is controlled by the geometrical and physical properties of the discontinuities. Numerical tools are commonly used for analysis of the stability of rock structures or evaluation of the failure process of rock masses. In this section, two models are specially designed and analyzed using the proposed approach to show the potential of the proposed approach in application of rock mechanics problems.

Frictional collapse of a rock mass system
One common application scenario of the numerical tools is the evaluation of collapse process of rock masses, such as the runout of rock avalanche and the impact of landslide to downstream structures or rivers. In this example, the proposed approach is used to investigate the influence of frictional properties on the collapse process of a rock block system. A jointed rock mass model is established as shown in Fig. 19. The region has dimensions of 20 m Â 50 m. Two sets of persistent joints cut this region into a rock block system. Dips and spacing of the two joints sets are 80 and 40 , 1:5 m and 2:5 m. It can be predicted that this block system will collapse with combination of joint geometry and low joint friction angle. The friction angle between rock blocks and boundary blocks is set to 35 . Two values are used for the friction angle between rock blocks: case 1: 15 ; case 2: 30 . Other physical parameters are as follows: block density 2700 kg/m 3 ; elastic modulus 2 GPa; Poisson's ratio 0:3; normal contact penalty 2 GPa; tangential contact penalty 1 Â 10 9 N/m; time step size 0:0001 s; and total simulation time 7:5 s. Comparison of the simulation results at 2:5 s, 5:0 s, and 7:5 s using different joint friction angles (case 1 of 15 and case 2 of 30 ) is shown in Fig. 20. A faster collapsing process and a larger runout mass were observed with lower friction angle. This example shows the proposed approach has potential in

Rock toppling and rock slumping
Another common application of numerical tools is to predict the potential kinetic failure pattern of rock slopes with certain physical parameters and geometrical patterns. In this example, the previous example [36] that involves both toppling and slumping simultaneously is used to show the accuracy of proposed approach in predicting the slope failure patterns. The model is shown in Fig. 21. The rock slopes in a valley region are cut by a joint set with dip angle of 75 and spacing of 3 m, and there is a sliding base with an inclined angle of 30 . The friction angle between the sliding base and the slope block is set to 32 , indicating that both the left slope and right slope will keep stable as one unique block. However, influenced by the joint set with a friction angle of 22 , the stable condition of both the left part and the right part will be broken. The left part will fail in slumping pattern, while the right part will fail with toppling pattern. Successful capture of simultaneous toppling and slumping behavior can verify the feasibility of the approach in modeling complex landslide behaviors. The above two-dimensional model and the extended threedimensional models have been used in verifications of discrete element approach [10] and discontinuous deformation analysis approach [63,65]. In this paper, the model is simulated with the proposed approach using the following parameters: block density 2700 kg/m 3 ; elastic modulus 2 GPa; Poisson's ratio 0:3; normal contact penalty 0:5 GPa; tangential contact penalty 1 Â 10 8 N/m; time step size 0:0001 s; total simulation time 3:5 s. The result after 3.5 s simulation is shown in Fig. 22. The proposed approach correctly predicted the slumping and toppling failure pattern. This shows the potential of the proposed approach in analysis of the stability analysis of rock block systems.

Conclusion
In this paper, the penalty function method with a shortestpath-based distance potential definition is proposed for contact force computation of two arbitrarily shaped convex polygons. This novel approach extended the definition of potential function for convex polygons and exhibit good robustness in treating irregular convex polygonal shapes. Currently, the proposed contact force computation approach is implemented in DDA framework with explicit representation of contact force. The contribution of this paper to discontinuous deformation analysis of convex polygon system is summarized as follows: (1) The boundary of entrance block is used to judge convex polygon overlap and detailed algorithm is given. A simple block cutting algorithm is also proposed to construct the intersection polygon when overlap occurs. These algorithms facilitate accurate and robust geometry analysis in contact of convex polygons. (2) The novel potential function definition for polygon lays a foundation for robust contact solution of irregular polygons. Specially for convex polygons, a simple algorithm is given for the computation of contact potential of any points in it. An integral The accuracy and robustness of the proposed approach are verified by the benchmarking examples. Moreover, its potential in kinetic analysis of complex polygonal block systems is shown by the two application examples. These advantages will promote the application of the proposed approach in any discontinuous methods that solve the interaction of discrete bodies with explicit representation of contact force.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.