Abstract
In the original discontinuous deformation analysis (DDA) method, the complete first-order displacement function is used to describe block movement and deformation, which induce constant stress and strain throughout the block. To achieve a more detailed stress distribution, Wachspress interpolation displacement function is employed to express the displacement of blocks in DDA, and the interactions between blocks are still governed under the original DDA. Displacements of the vertexes of all blocks constitute new freedom vectors, and the stiffness and force matrix formulations are derived again. In the new formulation, Wachspress interpolation ensures that the edges of the blocks are straight; therefore, contact detection can be processed based on the original DDA. Several classical examples are analyzed. The results show that the new formulation obtains similar configurations as the original DDA but provides more detailed and continuous stress distributions within block element.
1. Introduction
The discontinuous deformation analysis (DDA) method, which was proposed by Shi [1] in 1988, represents a novel method of analyzing the dynamic mechanical behavior of block systems under large displacement. As a theoretically rigorous numerical method, the DDA has been a focus of the discontinuous computing field since it was first proposed, and it has been used to address diverse problems in geotechnical engineering, such as landslide simulations [2], rock slope stability assessments [3], rockfall path tracking [4], blasting effect evaluation [5], stability analysis of tunnels [6], analysis of the shear creep deformation of structural planes [7], simulation of the fracturing process initiated by hydraulic pressure [8], crack propagation modeling [9], and seismic wave propagation modeling [10].
After numerous enhancements, the DDA has matured to solve more complicated problems. Solutions have been proposed to restrain false volume expansion when simulating large rotation, such as adopting the second-order approximation of and in the shape function [11], executing a postadjustment once the open-close iteration converges [12], introducing a new displacement variable based on the trigonometric function transformations [13] and accumulating the total strain components in fixed local frames [14]. To avoid the introduction of virtual springs, which are commonly used in classical DDA, the DDA has been successfully reconfigured using the Lagrange multiplier method [15], the complementarity theory [16], and the variational inequality theory [17]. Recently, a dual-form DDA has been proposed by Zheng [18], in which the contact forces instead of the block displacements are utilized as the basic variables and the compatibility iteration for the quasi-variational inequality guarantees the convergence and solution efficiency. To improve the efficiency of disk-based DDA, Beyabanaki and Bagtzoglou [19] presented a new contact model for nonrigid disks, in which disk-disk and disk-boundary contacts are transformed into point-to-line contacts. To simulate progressive failure, Jiang et al. [20] introduced a viscous damping component to absorb the kinetic energy of discrete blocks and defined convergence criteria for DDA solutions, which provided more objective standards for the application of the DDA.
Considering that 3D analyses are always preferred for practical engineering problems, many researchers are working on expanding the original 2D DDA to a 3D version. When the blocks have been identified via 3D discontinuities cutting [21], detecting contact rapidly and reliably appears to be the major barrier to developing an efficient 3D DDA. Contact algorithms, such as the incision body method [22], fast common plane method [23], and main plane method [24], have already been suggested. Recently, Shi [25] suggested a new contact theory named the E(A,B) algorithm, and it provides a more robust method of eliminating the greatest obstacle in 3D DDA programming.
In the original DDA, the first-order displacement function is adopted to describe the block movement and deformation, which induce constant stress and strain throughout the block. More exact stress distribution within block element is beneficial to the analysis of practical problems, and several approaches for determining such distributions have been attempted. One method is to assemble small blocks to form a larger block by adopting artificial joints and subblocks [26, 27]. When the subblocking method is used, the constraints between subblocks must be accurately imposed for reasonably reflecting the deformation of an actual block, which is difficult in practical applications. Another approach is to use high-order displacement functions in the DDA [1]. The second-order and third-order displacement functions have been successfully implemented [28, 29], and the ability to model complicated stress and strain fields have been verified. If high-order displacement functions are used, then the edges of a block are likely to be transformed to curves. Thus, the dissection of edges should be executed to exactly impose contact conditions, which is an additional troublesome job. In addition, scholars have added finite element meshes in blocks so that stress variations within the blocks can be evaluated [30, 31]. Considering that the finite element method (FEM) is familiar to researchers, this process is well understood. However, when the deformation and displacement of one block are illustrated by inner finite elements [32], the stress and strain fields are usually discontinuous within the block, which appears to represent a drawback compared with the adoption of high-order displacement functions in the DDA. To some degree, these attempts have successfully acquired more complicated stress and strain fields than the original DDA; however, extra tasks are likely involved or advantages of the original DDA may be lost.
Over the past two decades, the polygonal finite element method (PFEM) has achieved remarkable progress in theory and application [33, 34], and it provides new techniques for attaining a more detailed stress distribution of blocks in DDA. Some interpolations, e.g., Wachspress interpolation [35], Laplace interpolation [36], Sibson interpolation [37], and Maximum entropy interpolation [38], were proposed for constructing displacement functions in the PFEM. All displacement functions based on these interpolations have the ability to produce straight edges on polygons undergoing deformation. To begin with, Wachspress interpolation that originated from projective geometry defines the shape function at point using areas of triangular composed of point and the vertices. The gradient of the shape function can be derived directly, which facilities the calculation of strain and stiffness matrices for elastic analysis. Because the shape function is defined in the Cartesian coordinate system, the numerical integration could be performed on the physical elements. Then, both Sibson and Laplace interpolation are natural neighbor-based interpolations, which define the shape function at point using the information of Voronoi cells or Voronoi edges. For a better implementation, the isoparametric mapping and reference elements are usually involved in the computation [39]. Furthermore, enlightened by the Shannon entropy problem [40] in information theory, maximum entropy interpolation proposed by Sukumar [38] derives the shape function at point from a constrained optimization problem. Calculating the shape function requires solving an optimization problem by numerical iteration, which means that the subsequent computation is more complicated. By comparison, DDA users will have an easier access to Wachspress interpolation for little prior knowledge is involved, and Wachspress interpolation has a better performance being enrolled in the DDA codes. Therefore, Wachspress interpolation function is adopted in this paper to govern the displacement of each block in DDA. After contact detection in Cartesian coordinate system, the contact constraints between blocks are enforced as the original DDA. The displacements of vertexes of blocks constitute new freedom vectors and formulations for the modified DDA are derived. The C++ computer codes are developed for the new formulation and are used to solve several examples. The ability of the proposed scheme to obtain more detailed stress distributions within block elements is well verified by the calculated results.
2. Displacement Function of Blocks Based on Wachspress Interpolation
In the DDA method, the deformations and large displacements represent the accumulation of incremental displacements and deformations over many time steps. The movement and deformation of the th block within one time step are originally defined by six independent variables as follows:where and are the incremental displacement components of the centroid point of the th block in the horizontal and vertical directions, respectively; is the incremental rigid rotation angle around point ; and represents the incremental strains based on a small deformation assumption. The incremental displacement at any point in the th block can be related to : where is called the displacement transformation matrix of the block.Expression (2) has been verified to be the complete first-order displacement function [1]; therefore, each block presents constant stress and strain, which creates an obstacle to acquiring a more detailed stress distribution. Because most of the blocks in the DDA are polygons, Wachspress interpolation has great potential to be used developing a new displacement function for blocks. When Wachspress interpolation is employed to express the displacement of one block, the incremental displacements and of the vertexes form a new freedom vector of the block. Assuming that the th block contains vertexes within the current time step, could be rewritten as follows: where and are the incremental displacement components of the th vertex in the horizontal and vertical directions.
Wachspress interpolation requires that the shape function satisfies certain basic principles. For example, the five vertexes of a pentagon are arranged counterclockwise (shown in Figure 1). For the vertex , the edges including are called adjacent edges of , such as and , and the other edges are called opposite edges of , such as , , and .
The shape function must satisfy the following criteria:
(1) Continuous within the whole polygon;
(2) A function at vertexes, i.e.,
(3) Linear on adjacent edges
(4) Zero on opposite edges
Wachspress pointed out that a rational function could be adopted to build the required shape function and proposed an algebraic geometry approach to accomplishing this task [35]. Because the algebraic geometry approach is complex and impractical when the polygon has a large number of edges, some enhancements have been made to perfect the construction method of a shape function. The construction scheme proposed by Malsch and Dasgupta [39] is acknowledged as the simplest and most practical. Malsch and Dasgupta proposed that the shape function could be written as follows:where , , …, are the vertexes of a polygon arranged counterclockwise, is any point within the polygonal element, and A(·, ·, ·) is the area of a triangle determined by three points. An analysis of (6) indicates that the four basic principles are completely satisfied.
The incremental displacement at any point in the block can be expressed as follows:where is the shape function (6) associated with the th vertex. For a point on the boundary of a block, its displacement is actually a linear interpolation between the displacements of two relevant vertexes [35]. Adopting the displacement function (7), the edges of a block will still be straight lines after deformation, which is beneficial to contact detection.
Obviously, the shape function (6) can be replaced by a simpler formulation in the case of limited vertexes, such as a triangular finite element (n=3) or a quadratic finite element (n=4) with shape functions familiar to scholars. In this study, the shape function (6) is employed regardless of the number of vertexes because it has been verified that the element with shape function (6) will degenerate to those two classical finite elements if n=3 or n=4, which is beneficial to uniformly program codes.
Once the incremental displacement is obtained through the open-close iteration, the incremental strain vector within the current time step can be calculated as follows:in which is the strain matrix similar to the FEM and represents the relationship between strain and freedom. The incremental stress vector of the block is calculated byin which and is the elastic matrix of the th block. In the case of plane strain, iswhere is Young’s modulus and is Poisson’s ratio.
The total stress vector is accumulated as where represents the total stress vector of the block at the end of the mth time step and will be treated as the initial stress for the next time step.
3. Equilibrium Equations of the Block System
Because the DDA method applies the minimum total potential energy principle, for each block, the total potential energy is the summation of all potential energy sources, which include the elastic deformation of the block, initial stresses, point loads on the block, volume forces, and inertial forces. If one block contacts other adjacent blocks, the potential energy induced by contact forces or friction forces contributes to the total potential energy of the block.
Fixed points are abandoned in the new procedure. In the original DDA, fixed points are employed to enforce the displacement boundary conditions, and the potential energy induced by those constraint displacement points must be added to the total potential energy. However, because the new procedure considers the displacement of block vertexes as freedom vectors, displacement boundary conditions can be compulsively enforced by modifying the ultimate global stiffness matrix and load vector similar to the FEM; therefore, the potential energy associated with displacement boundary conditions is excluded from the total potential energy.
For a block system with blocks, to minimize the total potential energy, the equilibrium equations within the current time step can be expressed in the following matrix form:where is the displacement variables of block and is the generalized load vector acting on block , and these variables are represented by matrices ( represents the number of vertexes on block ). The stiffness matrix is a matrix depending on the material properties of block , and is a matrix defined by the contacts between block and .
In this section, based on the minimum total potential energy principle, the submatrices induced by all potential energy sources except the contact constraints are derived.
3.1. Submatrix of Block Elastic Deformation
The strain energy of the elastic stresses of the th block is derived as follows:In (14), the integration is performed over the entire area of the Ith block. For each time step, the blocks are assumed to be linearly elastic; therefore, the relationship between stress and strain can be denoted as follows:Therefore, the strain energy of (14) can be expressed as follows:If the strain energy is minimized by derivation, the stiffness matrix of the th block will be obtained and can then be added to the global stiffness matrix:
3.2. Submatrix of Initial Stress
The accumulated stresses over the previous time steps will be transferred to the next step as the initial stress loading, which means that the stress in (12) will be considered the initial stress for time step m+1. The potential energy of the initial stress is given as follows:If is minimized by derivation, the following vector will be obtained and can then be added to the vector in the global force vector:
3.3. Submatrix of Point Loading
Assuming that the point loading force acts on any point within the th block, the potential energy of this point loading is as follows:If is minimized by derivation, then the following vector will be obtained and can be added to the vector in the global force vector:
3.4. Submatrix of Volume Forces
Assuming that represents the constant body forces acting on the entire volume of the th block, the potential energy induced by body forces is as follows:If is minimized by derivation, the following vector will be obtained and can be added to the vector in the global force vector:
3.5. Submatrix of Inertial Forces
In the DDA, the time steps are used in both static and dynamic analyses. The static analysis assumes that the initial velocity at the beginning of the current time step is 0. When the dynamic analysis is executed, the velocity at the end of the previous time step is used as the initial velocity at the beginning of the current time step. The term (), )) denotes the time-dependent displacement of any point of the th block and denotes the mass per unit area; thus, the inertial force acting on a per unit area is as follows:The potential energy of the inertial force is given as follows:The Newmark scheme is adopted to dispose the accelerated velocity as the original DDA. We haveandwhere is the time interval of the present time step and is the velocity at the end of the previous time step.
3.6. Numerical Integration Scheme
Collecting the contributions of all potential energy sources mentioned above, we haveIn (28) and (29), the components and are actually rational functions, which means that the integration cannot be accomplished analytically. To overcome this problem, the original polygonal block is decomposed to several triangles, and a Hammer integration scheme is employed. The details of the integration scheme are given as follows.
As shown in Figure 2(a), for any given polygon, the integration on the whole area can be obtained as follows:For simplicity, each subtriangle is formed by two adjacent vertexes and the centroid point of the original polygon. The Hammer integration provides a simple method of calculating the integration of each subtriangle. According to the Hammer integration, where , , and are the area coordinates of an integration point in the triangle and is the number of integration points within one triangle. If m=4, the numerical integration precision of the Hammer integration will reach . Four integration points are employed as in Figure 2(b), and their area coordinates and related weighting coefficients are listed as follows:where S represents the area of the triangle.
(a) Dividing a polygon into small triangles
(b) Four Hammer points in a subtriangle
4. Submatrices of the Contacts between Blocks
The advantage of the DDA is the rigorous contact model that governs the interactions of blocks, and the model includes two main steps: “contact detection” and “contact mechanics”. When the displacements of blocks are governed by (7), the edges of blocks remain straight from beginning to end; hence, “contact detection” can be processed as the original DDA. However, because the freedom vectors have been changed in the new formulation, the “contact mechanics” should be adjusted accordingly.
No interpenetration and no tension are required on the contact surface when two adjacent blocks contact. To achieve the contact conditions, stiff springs are applied along normal and tangential directions. In this section, the submatrices related to contact are derived again.
4.1. Normal Stiffness Submatrix
As shown in Figure 3, a vertex and an entrance line constitute a classical contact pair in DDA. The coordinates of point , , and are denoted as , , and , respectively, and the displacements are , , and in the current time step, respectively. Considering that displacements in each time step are required to be a small value, the first-order infinitesimal value is retained and the second-order infinitesimal value is neglected. The distance from vertex to entrance line can be expressed as follows:where andFurthermore, is expressed as a formulation about , , and : in whichTo satisfy the no interpenetration requirement, the normal spring is added when d<0, and the deformation energy of the normal spring is as follows: where is the stiffness of the normal spring. By minimizing the deformation energy, nine 2×2 submatrices are obtained and can be added to the global stiffness matrix as shown in (39) and three 2×1 submatrices are obtained to be added to the global load vector as shown in (40):where is the generalized load vector acting on vertex . This vector is a 2×1 vector. The stiffness matrix is a 2×2 matrix related to vertex , and is a 2×2 matrix caused by contact restraints.
4.2. Tangential Stiffness Submatrix
As shown in Figure 4, the same contact pair is analyzed in the tangential direction. The point is the lock point of vertex at the entrance line . The coordinate and displacement of are denoted as and . Similarly, by neglecting the second-order infinitesimal value, we generate the relative tangential distance along line :whereand is the same as in (35).
As mentioned in the second section, the displacement of point can be obtained by a linear interpolation between and : Therefore, is illustrated as a formulation about , , and :in whichWhen the contact pair is locked in the tangential direction according to Mohr-Coulomb friction law, the tangential spring is added and the related deformation energy is as follows: where is the stiffness of the tangential spring. By minimizing the deformation energy, nine 2×2 submatrices are obtained and can be added to the global stiffness matrix as shown in (47) and three 2×1 submatrices are obtained and can be added to the global load vector as shown in (48):
4.3. Friction Force Submatrix
As shown in Figure 5, the contact pair is reanalyzed in the tangential direction. When the vertex slides along the entrance line , the tangential spring will be replaced by a pair of friction forces. The direction of friction is defined by the relative tangential displacement between the master point and its lock point . Assuming that is the normal spring stiffness and is the normal interpenetration, the friction force should be where is the friction angle and is a sign function about : Figure 5 provides a rule to determine whether x is positive. With this rule, the friction direction can be denoted as follows:When the value and direction of the friction force are ascertained, they can be treated as a pair of point loads acting at point and point . With the new displacement function, for one point located on the boundary, a point load will result in the load vector of two associated vertexes.
Finally, the contribution of friction force can be expressed as three 2×1 submatrices that can be added to the global load vector as shown inThe global equations (13) are assembled by collecting the submatrices of all contacts between blocks and from (28) to (29) for all blocks, and they are further modified according to existing displacement boundary conditions. After solving the global equations, the states of each contact are inspected and the springs are adjusted to reflect the difference between the surmised contact states and the calculated contact states, which is called the “open-close” iteration. When the surmised contact state is consistent with the calculated contact state, the computation terminates and the freedom vectors are finally determined.
5. Numerical Examples
5.1. Sliding Problem
As shown in Figure 6, a rectangular block is laid on a ramp with a slope angle of 45°. The block and ramp have the same material properties: density ρ=2.5×103 kg/m3, Young’s modulus E=35 MPa, Poisson’s ratio ν=0.3, and gravity acceleration g=9.8 m/s2. The external load is only generated via gravitational force. The three vertexes of the ramp are all fixed as the displacement boundary conditions.
Assuming the friction angle on the interface is , the block will slide along the ramp when φ<45°; in addition, its velocity has an analytical solution:where is the sliding velocity (m/s) of the block and is the elapsed time (s).
Four friction angle cases are studied: 0°, 15°, 30°, and 40°. Let the time step length Δ=0.02s and the max step displacement ratio δ=0.10; subsequently, 35 time steps are calculated. Figure 7 shows a comparison between the analytical solutions and the results from the proposed scheme. The velocities from the new scheme are consistent with the analytical solutions for all cases of friction angles.
The original DDA will provide high-precision results similar to this example, and the original DDA and the new formulation induce different stress distributions within the block. For this example, the stress results obtained by the original DDA and the proposed scheme are listed in Table 1. When the friction angle is 0°, the results of the proposed scheme are basically consistent with the original DDA because in this case, the block has a constant stress distribution, which conforms to the assumption in the original DDA. However, the results of the proposed scheme are quite different from that of the original DDA in the case of friction force. An analysis of the data indicates the stress results of the original DDA closely approximate the average value of the stress results provided by the proposed scheme.
Figures 8 and 9 show the principal stress distributions of the model with the four friction angles. In the case of friction angle 0°, the principal stresses and are roughly constant in the block. Principal stress is positive and is negative, which indicate that the block is in tension along the direction of the principal stress and compressed along the direction of .
(a) Friction angle 0°
(b) Friction angle 15°
(c) Friction angle 30°
(d) Friction angle 40°
(a) Friction angle 0°
(b) Friction angle 15°
(c) Friction angle 30°
(d) Friction angle 40°
When the friction angle increases, the friction force acting along the bottom surface of the block gradually increases, which leads to apparent changes in the stress distributions. When friction force occurs, the minimum value of stress is negative and remains at the bottom right corner, which means that this area is compressed. With increases in the friction angle, the minimum value of stress decreases accordingly, and the compressive region spreads from the right to the left. When the friction angle is 40°, the regions close to the contact surface are almost compressed along the direction of principal stress . The maximum value of stress rises as the friction angle increases, which results in the distortion and deflection of the stress contour.
Principal stress is negative throughout the block in most cases, which represents the compressive state within the entire block. The minimum value of principal stress is located at the bottom right corner and falls rapidly as the friction angle increases. The maximum value of occurs at the top left corner, and the value increases when the friction angle increases. When the friction angle reaches 40°, the maximum value of is positive, which indicates that the top left corner is under tension.
5.2. Arch Bridge Problem
The arch bridge problem is a classical model designed by Shi and often used to test the validity of procedures based on the DDA. Figure 10 shows an arch bridge that consists of seven blocks. Each block has the same shape and area (0.17227 m2), and four vertexes located at the leftmost and rightmost boundary are fixed as displacement boundary conditions. A vertical downward sustaining point load P=0.05 kN acts on the central block. All blocks have the same property with the unit weight γ=2.0 kN/m3, Young’s modulus E=10 MPa, and Poisson’s ratio ν=0.2. Using the time step length Δ=0.01 s and the max step displacement ratio δ=0.03, we conduct a dynamic analysis of the example using the original DDA and the proposed formulation.
The configurations and stress distributions of the arch bridge for the original DDA and the new formulation after 30 time steps are shown in Figure 11, which assumes that no friction occurs between blocks. The results show that the new formulation generates a similar configuration as the original DDA. However, stress distributions from the new formulation are more explicit than those of the original DDA.
(a) Stress from the original DDA
(b) Stress from the original DDA
(c) Stress from the new formulation
(d) Stress from the new formulation
The stress distributions of the blocks fall into two categories. First, the stresses have significant variations in the leftmost and rightmost blocks. The values of stress and are both negative in the regions close to fixed boundaries and change from negative to positive when the measured points move toward the center of the arch. For the boundary blocks, the regions close to fixed boundaries suffer biaxial compression, whereas the regions adjacent to the middle blocks are actually in tension. Second, the variations of all stresses are small in the middle blocks. An analysis of the stress distribution indicates that only the stress state of the central block is consistent with the original DDA. Along the direction of stress , the bottom regions are in compression in the other middle blocks, whereas along the direction of stress , some top regions are in tension. Finally, the results of the original DDA are still close to the average value of the new formulation.
Then, the above example is reanalyzed under the following assumptions: friction angle =30° and cohesion c=0. The arch bridge is stable after 30 time steps, and the vertical load acting on the central block is distributed to the adjacent blocks one by one. The configurations and stress distributions from the original DDA and the new formulation are illustrated in Figure 12.
(a) Stress from the original DDA
(b) Stress from the original DDA
(c) Stress from the new formulation
(d) Stress from the new formulation
The two methods both obtain a negative stress throughout the model, and its absolute value is larger than stress , which indicates extrusion between blocks. Compared with the original DDA, the new formulation reveals additional details about the extrusion effect. In the middle blocks, stress reaches its minimum value at the bottom; however, stress reaches its minimum value at the top in the leftmost and rightmost blocks. Thus, the extrusion effect varies with the position of a block in this arch model. Combined with stress , most of the blocks are not in biaxial compression everywhere. For a boundary block, stress changes from negative to positive when the measurement point moves toward the center. For the central block, stress changes from positive to negative when the measurement point moves downward. An analysis of the joint surfaces shows that stress always has a different sign on both sides.
The difference between the two methods is not remarkable in Figure 12, which is caused by the stress tag setting. For ease of comparison, all panels in Figures 11 and 12 adopt the same stress tag; however, the stress distribution obtained by the new formulation is more complicated. Figure 13 shows the stress distribution with an independent tag.
(a) Stress
(b) Stress
5.3. Voronoi Polygons Problem
Voronoi polygonal structures are commonly used to simulate the microstructures of crystals and other materials, and they present considerable advantages for creating complex particle shapes and sizes. Generally, the particles are typical polygons, and for each particle, the number of edges may be different. The new formulation is preferable for simulating block models generated by the Voronoi map method.
As shown in Figure 14, the block model with a Voronoi polygonal structure is laid on a ramp. The ramp and blocks have the same properties, with a unit weight γ=21 kN/m3, Young’s modulus E=10 GPa, and Poisson’s ratio ν=0.25. The mechanical parameters of the joints throughout the model have a friction angle of φ=40° and cohesion of c=1 MPa. The contact surface between the blocks and the ramp uses the same parameters applied to the joints. Four vertexes of the ramp are fixed as displacement boundary conditions.
Assuming that the external load is only generated via gravitational force, we conduct a static analysis of the example with the original DDA and the new formulation, and it includes a time step length of Δ=0.01 s and a max step displacement ratio of δ=0.03. The configurations and stress distributions after 55 time steps are illustrated in Figure 15.
(a) Stress σ1 from the original DDA
(b) Stress σ2 from the original DDA
(c) Stress σ1 from the new formulation
(d) Stress σ2 from the new formulation
The two methods obtain a negative stress throughout the model, and the absolute value increases when the measurement point moves downward because only gravity guides the stress distribution and deformation. Compared with stress , stress has small absolute values and variations in the model. Except for some regions close to joints, the results from the new formulation reveal that stress ranges from -18.54 N/mm2 to 18.04 N/mm2.
Compared with the original DDA, the new formulation shows issues with the stress concentration around some vertexes. In the original DDA, the contact forces are developed by mutual movements between blocks regardless of whether a static or dynamic analysis is performed. This approach is not changed in the proposed formulation, which results in concentrated forces acting on vertexes. Therefore, the stress concentration issue is reasonable.
A horizontal rightward sustaining point load P=200 kN is added to act on the upper left block in Figure 16. The example is reanalyzed with the same calculation parameters and time steps as before.
Figure 17 shows the configurations and stress distributions from the original DDA and the new formulation. The distributions of stress and stress have clearly changed, especially for the new formulation. Most regions are in biaxial compression under a combination of a horizontal pushing force and gravitational force. The upper left block presents the maximum pressure stress . Because of the effect of the horizontal pushing force, the regions with a positive stress changed and the area decreased. Stress increases around some joints because of the transmission of the horizontal pushing force.
(a) Stress from the original DDA
(b) Stress from the original DDA
(c) Stress from the new formulation
(d) Stress from the new formulation
6. Discussion
The new formulation improves upon existing approaches because the complicated and continuous stress distributions within block elements are provided and all edges of blocks remain straight. Regardless of whether artificial joints or finite element meshes are involved in the block, the block dissection should be executed and the stresses will be discontinuous within the block. The new formulation avoids block dissection and generates a continuous stress distribution. When higher-order displacement functions are adopted in the DDA, the edges of blocks are always transformed from lines to curves, which induce problems for contact detection. To satisfy the contact constraints between blocks, the edges should be dissected and then contact detection should be carefully performed. In the new formulation, the edges of blocks remain straight lines, which means that contact detection can be processed according to the original DDA. Of course, the DDA with higher-order displacement functions presents certain remarkable abilities; e.g., one straight line polygon can be transformed into one curved polygon. However, for most rocks, the strain in the block is small, even if the block is on the point of rupture. Accordingly, the deformation of a block is small compared with the movement, and the assumption that the edges of a block are still straight is acceptable.
The limitation of the new formulation is that only those blocks with the shape of a convex polygon can be included in the block system because the displacement function based on the Wachspress interpolation is not suitable for concave polygons [33]. The displacements of block vertexes form new freedom vectors in the proposed formulation, which implies that, for one block, the detailed level of stress distribution relies on the number of its vertexes rather than the number of Hammer points involved in the integration. The extreme case is a block with three vertexes, for which only constant stresses can be gained by the proposed formulation. Because the Wachspress interpolation scheme is effective for all nonconcave polygons, the issue can be solved by appending extra vertexes on edges. However, introducing additional vertexes to generate a detailed stress distribution is not feasible in the new formulation because the m+1th time step is always calculated based on the configuration results at the end of the mth time step. If one additional vertex is added on one edge, then the original nonconcave polygon may be transformed to a polygon containing a concave angle after deformation and the computation will fail. A more robust displacement interpolation scheme is expected to be developed for the uniform processing of various polygons, which will remedy the limitations of the proposed formulation.
7. Conclusions
To acquire a more detailed stress distribution, the displacement function of block in DDA is modified in this paper. For each block, the displacements of its vertexes constitute new freedom vectors, the displacement is expressed based on Wachspress interpolation function, and formulations of stiffness and force matrices are then derived accordingly. The contacts between blocks are still governed as in the original DDA; thus, contact detection can be processed based on the original DDA because the edges of blocks are always straight. Three examples are analyzed using the original DDA and the proposed formulation. The results show that the new scheme has similar configurations as the original DDA but generates more detailed and continuous stress distributions within block elements. The new formulation has greater advantages because it provides complicated and continuous stress distributions while avoiding extra work in the contact analysis. Because of the limitations of existing Wachspress interpolation approaches, the new formulation is not suitable for concave polygons and lacks the ability to obtain detailed stress distributions of blocks with limited vertexes. An advanced displacement interpolation scheme suitable for uniformly processing various polygons is expected, and it will improve the robustness of the new formulation.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study is supported by the Open Research Programme of the Hubei Key Laboratory of Disaster Prevention and Mitigation (under Grant no. 2017KJZ03) and National Natural Science Funds (Project no. 51409150).