Vertex Displacement-Based Discontinuous Deformation Analysis Using Virtual Element Method

To avoid disadvantages caused by rotational degrees of freedom in the original Discontinuous Deformation Analysis (DDA), a new block displacement mode is defined within a time step, where displacements of all the block vertices are taken as the degrees of freedom. An individual virtual element space V1(Ω) is defined for a block to illustrate displacement of the block using the Virtual Element Method (VEM). Based on VEM theory, the total potential energy of the block system in DDA is formulated and minimized to obtain the global equilibrium equations. At the end of a time step, the vertex coordinates are updated by adding their incremental displacement to their previous coordinates. In the new method, no explicit expression for the displacement u is required, and all numerical integrations can be easily computed. Four numerical examples originally designed by Shi are analyzed, verifying the effectiveness and precision of the proposed method.


Introduction
Discontinuous Deformation Analysis (DDA) [1], a novel numerical method for analyzing the dynamic mechanical behavior of a block system in cases of large displacement, was verified as an effective tool in solving a variety of discontinuities in rock problems [2,3]. In rock-mass engineering, DDA has been employed to handle a great deal of problems, e.g., landslide process simulations [4,5], slope stability assessments [6,7], blasting effect evaluation [8], crack propagation simulation [9,10], seismic wave propagation analysis [11], and rock burst prediction [12]. Applications of 2D DDA in the modeling of rock-mass dynamics until 2017 were summarized by Ning [13].
Numerous enhancements have been put forward to improve the performance of traditional DDA. To alleviate the sensibility of penalty parameters, contacts between the blocks were remodeled using the Lagrange multiplier method [14], Augmented Lagrangian method [15], Complementary theory [16], and Variational Inequality theory [17]. Apart from modifying kinetic velocity using the dynamic factor [18], the damping effect was used to reflect the energy dissipation by imposing viscous boundary conditions [19] and by adding viscous forces into the equilibrium equation [20]. A couple of strategies, such as adopting a higher-order displacement function [21] and partitioning blocks using finite element mesh [22], were introduced to acquire a more detailed stress field in each block. Moreover, many efforts were made to develop a 3D version of DDA [23][24][25][26], considering that 3D DDA is preferred for problems in practical engineering. The biggest bottleneck to establishing a robust 3D version is the lack of an excellent contact theory [27]. Shi [28] proposed the entrance block theory to facilitate contact treatment in 3D DDA. In addition, some efforts, such as GPU-based parallel computation [29] and explicit computation [30], were made to improve the computational efficiency when simulating large-scale problems.
Within a time step, the degrees of freedom d in the traditional DDA were defined for a block by six independent variables, which consists of two incremental rigid translations, one incremental rigid rotation angle, and three incremental constant strain components. For a point x in the block, the incremental displacement u is calculated as the product of d and the transformation matrix T; then, the global equilibrium equation is derived from this displacement function. The displacement function uses the first-order approximation of sin r 0 = r 0 and cos r 0 = 1 for a small rotation angle r 0 , which induces accumulated errors over the time steps and sometimes causes issues. The most noticeable issue is the unreasonable volume expansion when a large rotation occurs. To restrain this issue, some scholars proposed various emendations [31][32][33], in which the post-adjustment strategy was most popular in the existing DDA codes [34,35]. This method, albeit simple and effective, resulted in a displacement different from the displacement evaluated in the equilibrium equation. Even with a tiny difference, the contact state might be entirely changed for a contact pair. Then, to simulate the continuous-discontinuous deformation, numerical methods for the continuum and non-continuum models were increasingly coupled. DDA was coupled with finite element method (FEM) by some scholars [9,22] for simulation of a rock failure process. In classical FEM, displacements at the element nodes are chosen as the degrees of freedom. This inconsistency in the degrees of freedom between DDA and FEM causes a barrier in developing a compact and efficient code for the coupling method.
In a recent study [36], the displacement of blocks in DDA was reformulated by using Wachspress interpolation to achieve a higher-order stress distribution within blocks. This approach selected displacements at the vertices as the degrees of freedom for a block. In this study, a DDA with such new degrees of freedom is called a vertex displacement-based DDA. This method was improved using the Polygonal Finite Element method (PFEM). Besides Wachspress interpolation [37], a series of interpolations [38][39][40] were proposed to construct the displacement functions for a polygonal element with more than four nodes in PFEM. By taking displacements at the vertices as the degrees of freedom, these displacement models can avoid demerits caused by the degrees of freedom in an original DDA. However, because the rational function is involved in the shape function, the theoretical derivations and calculations of the stiffness and mass matrix are quite complicated and it is hard to ensure precision of the numerical integration [41,42].
In 2013, the Virtual Element Method (VEM) was proposed to handle the very general polygons, dispensing sophisticated integrations on the element [43][44][45][46]. VEM was considered the evolution of the mimetic finite difference method and a generalization of FEM. For flexibility with regard to mesh generation and element shapes, VEM has become a hot topic in numerical methods since it was proposed [47][48][49][50][51][52]. In this study, a new vertex displacement-based DDA is developed using the strength of VEM. Defining the degrees of freedom by using the incremental displacements u at the vertices of a block, an individual virtual element space V 1 (Ω) is adopted to describe the displacements of points in the block, and the projector Π P u from V 1 (Ω) on the linear displacement space P 1 (Ω) is deduced. Next, the total potential energy is investigated for the block system. In the potential energy, the bilinear forms of u are expressed as the summation of the exact solution of Π P u and an approximation of u-Π P u. The potential energy induced by the contact restraints are derived using the new degrees of freedom. Then, for the block system, the global equilibrium equation is derived based on the principle of minimum potential energy. Finally, the open-close iteration strategy is employed to resolve the global equilibrium equation as the original DDA. The proposed method avoids the issues attributable to first-order approximation for a small rotation angle r 0 in the original DDA and has a higher computational efficiency than vertex displacement-based DDA using displacement functions in PFEM. The validity and effectiveness of the proposed method are verified by several numerical examples.

Basic Principles of DDA
As Figure 1 shows, a DDA block system always consists of n b discrete blocks with their individual domains and boundaries. The domain Ω I of block I is bounded by ∂Ω I , which is usually composed of the Dirichlet boundary Γ I u and the traction boundary Γ I t . The displacement on Γ I u is prescribed as û I , and the surface traction on Γ I t is denoted by t I . Here, the superscript I is the block index. The deformations and large displacements of a block in DDA are accumulated by the incremental displacements and deformations over time steps. The displacements of blocks are independent from each other, and contact constraints are imposed on the interactions between blocks.
angle r0 in the original DDA and has a higher computational efficiency than vertex displacement-based DDA using displacement functions in PFEM. The validity and effectiveness of the proposed method are verified by several numerical examples.

Basic Principles of DDA
As Figure 1 shows, a DDA block system always consists of nb discrete blocks with their individual domains and boundaries. The domain Ω I of block I is bounded by ∂Ω I which is usually composed of the Dirichlet boundary Γ I u and the traction boundary Γ I t The displacement onΓ I u is prescribed as û I , and the surface traction on Γ I t is denoted by t I Here, the superscript I is the block index. The deformations and large displacements of a block in DDA are accumulated by the incremental displacements and deformations over time steps. The displacements of blocks are independent from each other, and contac constraints are imposed on the interactions between blocks. Within a time step, DDA requires the incremental displacements u = (u 1 , u 2 , …, u nb to be the unique minimizer of the total potential energy of the block system: where J'(v) represents the energy due to the contact constraints. A key feature of DDA is that rigorous contact constraints are used to manage the interactions between blocks. For the contact pair marked in Figure because all blocks are physically isolated in DDA. For ease of description, the superscrip indicating blocks is omitted unless otherwise noted. The bilinear form α (u, v) for u∈V and v∈V represents the energy due to elastic deformation and the inertial force: where ε is the incremental constant strain decided by u: Within a time step, DDA requires the incremental displacements u = (u 1 , u 2 , . . . , u nb ) to be the unique minimizer of the total potential energy of the block system: in which V is the incremental displacement space of the block system: where J'(v) represents the energy due to the contact constraints. A key feature of DDA is that rigorous contact constraints are used to manage the interactions between blocks. For the contact pair marked in Figure 1, J'(v) induced by the contact constraints is determined by u I ∈V I and u L ∈V L . The bilinear form α(u, v) and the linear form f (v) are computed for u∈V and v∈V by because all blocks are physically isolated in DDA. For ease of description, the superscript indicating blocks is omitted unless otherwise noted. The bilinear form α (u, v) for u∈V and v∈V represents the energy due to elastic deformation and the inertial force: where ε is the incremental constant strain decided by u: in which ∇ is the gradient operator. σ is the incremental Cauchy stress related to ε by the constitutive equation: where D is a constant elasticity tensor. ρ is the mass per unit area. For dynamic analysis, the inertial force is essential in DDA. New-mark time integration scheme is adopted in the original DDA, and the acceleration is assumed to be constant within a time step. Denoting by V 0 the velocity of block I at the origin of a time step and by ∆, the time interval of the time step, the acceleration is computed as The linear form f (v) is defined by where b denotes the constant body force and σ 0 is the constant stress accumulated over the previous time steps. When the displacement function is determined for a block with the associated degrees of freedom, Equation (1) induces the global equilibrium equations with the following matrix formulation for the block system.
in which d I is the degrees of freedom concerning block I. F I is the generalized load vector with the same dimension of d I . Denoting the dimension of matrix denoting the contact restraints on blocks I and L. If no contact pair is provided by blocks I and L in the current time step, K IL is zero.

Demerits Caused by the Original Degrees of Freedom and Previous Attempts to
where u 0 and v 0 represent the increments of the rigid horizontal and vertical translations, respectively; r 0 denotes the increment of the rigid rotation angle around the central point (x 0 , y 0 ) of block Ω; and (ε x , ε y , γ xy ) denotes the increments of the constant strains. Under small deformation assumption, the increments of the displacement u = (u, v) T at a point x = (x, y) in block Ω are computed as u = Td (11) and T is the translation matrix as follows: Equation (11) is a variation of a standard first-order displacement function [1]. The displacement caused by rigid movement and the deformation of the block constitute the displacement of a block under small deformation assumption. Therefore, the displacement function for finite rotation should be Comparing Equation (13) with Equation (11), the approximations of cos r 0 = 1 and sin r 0 = r 0 are adopted for the small rotation angle r 0 in Equation (11). The approximation errors accumulated over time steps can lead to false volume expansion when a large rotation occurs. Various modifications were suggested to remedy this defect, in which the post-adjustment strategy is most popular in the existing DDA codes. After vector d is obtained, the post-adjustment strategy employs Equation (13) to calculate the incremental displacement u. Although simple and effective in most cases, the resulting displacement must be different from the displacement estimated in the equilibrium equation. Sometimes, for a contact pair, the state adopted in the equilibrium equation may be not coincident with the geometry relationship of the resulted configurations. An example is provided in Section 5.1 to demonstrate this issue.
Besides that, it is noticed that DDA was coupled with a range of numerical methods for continuum models to simulate the continuous-discontinuous deformation, e.g., crack propagation [9,22]. Most popular numerical methods in terms of continuum mechanics select the displacements at the nodes of an element as the degrees of freedom. DDA defines the degrees of freedom using the constant strain, the rigid translation, and the rigid rotation referring to the centroid of a block. The difference in the degrees of freedom causes a barrier in developing a compact code when coupling DDA to other numerical tools.

Previous Attempts to Construct Vertex Displacement-Based DDA
By taking the increments of the displacements at the vertices as the degrees of freedom, vertex displacement-based DDA can remedy the demerits induced by the original degrees of freedom. Supposing that a block Ω has n vertices, the incremental displacements u i = (u i , v i ) T at the vertex x i constitute the new degrees of freedom: In this way, the penetration value of a contact pair in the equilibrium equations is directly calculated within a time step and the coordinates of the vertices are updated by directly adding the incremental displacements to their previous coordinates, inducing penetrations in the equilibrium equation and allowing the updated configurations to coincide.
By using the new degrees of freedom, the increment u of the displacement at a point x in a block Ω can be computed in a similar method to FEM, that is (15) where N i is the shape function in terms of the vertex x i . Because the block in DDA might have more than four vertices, the shape functions developed in PFEM were introduced to construct N i [36]. In the recent study [36], Wachspress interpolation was selected as the shape function for developing vertex displacement-based DDA. Sukumar and MalsCh concluded the properties required for N i to develop an effective PFEM [38]. As shown in Figure 2, a block Ω has n vertices located at x i and its boundary is composed of n straight edges. The red edges with vertex x i are denoted by Γ i . N i should be satisfied.
(1) Partition of unity, boundedness, and nonnegativity: (2) Kronecker-delta property: (4) With C ∞ being within block Ω while C 0 is on the boundary, N i must be piece-wise linear along Γ i but vanishes at the other edges. This property ensures that the boundary of a straight line is still a straight line after deformation, which benefits contact detection and contact condition imposition.
(2) Kronecker-delta property: (17) (3) Linear completeness: (18) (4) With C ∞ being within block Ω while C 0 is on the boundary, Ni must be piece-wise linear along Γi but vanishes at the other edges. This property ensures that the boundary of a straight line is still a straight line after deformation, which benefits contact detection and contact condition imposition. In the last two decades, generalized barycentric coordinates were successfully adopted as the shape function in PFEM. A series of generalized barycentric coordinates, such as Wachspress interpolation [37], maximum entropy interpolation [38], mean value coordinates [39] and Harmonic coordinates [40], were developed, marking remarkable progress in the theory and application of PFEM. Take Wachspress interpolation for example; the shape function Ni is defined for a point x within a polygon as follows. (19) where x1, x2, …, xn are the vertices of the polygon arranged counterclockwise and A(xi, xj, xk) represents the area of a triangle connecting three points, such as in Figure 3. By observation, the basic principles required for Ni are completely satisfied by Definition (19). However, the definitions of Ni in PFEM are the rational functions, which make the derivation and integration process more cumbersome when computing the block matrices. Though some enhancements [38] were made to simplify the formulas of Ni, the rational function is indispensable to definite Ni, which creates difficulties in the integration over a polygon. Up until now, the most popular way to solve this issue was partitioning the polygon into subdomains with certain shapes. Then, the Hammer or Gauss integration scheme was employed to compute the integration on each sub-triangle or sub-rectangle. The summation of the integrations on subdomains was considered the integration over the original polygon. While the integration error was considerable, some remedies were proposed to address the integration error for PFEM [42]. In conclusion, the additional work in the derivation and programming is very expensive when using In the last two decades, generalized barycentric coordinates were successfully adopted as the shape function in PFEM. A series of generalized barycentric coordinates, such as Wachspress interpolation [37], maximum entropy interpolation [38], mean value coordinates [39] and Harmonic coordinates [40], were developed, marking remarkable progress in the theory and application of PFEM. Take Wachspress interpolation for example; the shape function N i is defined for a point x within a polygon as follows.
where x 1 , x 2 , . . . , x n are the vertices of the polygon arranged counterclockwise and A(x i , x j , x k ) represents the area of a triangle connecting three points, such as in Figure 3. By observation, the basic principles required for N i are completely satisfied by Definition (19).  (15) is adopted to govern the block displacement, some advantages of the original DDA can compromised.
x i Figure 3. Triangles involved in the calculation of Ni using Wachspress interpolation.

Virtual Element Spaces
VEM [43][44][45][46] provides a new way to construct a vertex displacement-based DDA. Because DDA blocks are physically isolated, an individual virtual element space Vk(Ω) can be defined for a block as follows: in which Γ denotes a edge of Ω, k ≥ 1 is a fixed integer index indicating the order of accuracy of the approach, and Pk(Ω) is the space of polynomials of degree less than or equal to k in Ω. Obviously, [Pk(Ω)] 2 ⊆Vk(Ω). For simplicity, [Pk(Ω)] 2 is denoted as Pk(Ω). However, the definitions of N i in PFEM are the rational functions, which make the derivation and integration process more cumbersome when computing the block matrices. Though some enhancements [38] were made to simplify the formulas of N i , the rational function is indispensable to definite N i , which creates difficulties in the integration over a polygon. Up until now, the most popular way to solve this issue was partitioning the polygon into subdomains with certain shapes. Then, the Hammer or Gauss integration scheme was employed to compute the integration on each sub-triangle or sub-rectangle. The summation of the integrations on subdomains was considered the integration over the original polygon. While the integration error was considerable, some remedies were proposed to address the integration error for PFEM [42]. In conclusion, the additional work in the derivation and programming is very expensive when using the interpolation function N i in PFEM. When Equation (15) is adopted to govern the block displacement, some advantages of the original DDA can compromised.

Virtual Element Spaces
VEM [43][44][45][46] provides a new way to construct a vertex displacement-based DDA. Because DDA blocks are physically isolated, an individual virtual element space V k (Ω) can be defined for a block as follows: (20) in which Γ denotes a edge of Ω, k ≥ 1 is a fixed integer index indicating the order of accuracy of the approach, and P k (Ω) is the space of polynomials of degree less than or equal to k in Ω. Obviously, [P k (Ω)] 2 ⊆V k (Ω). For simplicity, [P k (Ω)] 2 is denoted as P k (Ω).
We take k = 1 in this study. There are three main reasons for this choice: (i) The degrees of freedom for V 1 (Ω) only involve the values of v at the vertices of Ω, in accordance with the new degrees of freedom defined in Equation (14). (ii) The incremental displacements are expected to be piece-wise linear at the block boundary, which is satisfied by V 1 (Ω). (iii) In the view that the original DDA adopted a complete first-order displacement function, the degree of accuracy is not threatened by using V 1 (Ω). Therefore, an individual virtual element space V 1 (Ω) is defined for a block, and the associated degrees of freedom are the new degrees of freedom defined in Equation (14). The dimension of the space V 1 (Ω) is 2n, where n denotes the number of vertices for the block.
Provided that there are n b blocks in the block system, the total virtual element space can be given as and the total degrees of freedom comprise the increments of displacements at the vertices of all blocks. As a result, the solution space V in Equation (1) is replaced by V 1 .

The Projection Operator
To minimize the total potential energy of the block system, the computation of α (u, v) from V 1 (Ω) × V 1 (Ω) to R defined in Equation (4) is inevitable, which can be accomplished in the framework of VEM theory. Referring to [47], two projectors Π ∇ u: V 1 (Ω)→P 1 (Ω) and Π 0 u: V 1 (Ω)→P 1 (Ω) are respectively defined using the orthogonality conditions: and where Π 0 u cannot be obtained from (23) because the zero-and first-order moments of u in Ω are unknown. To solve this issue, an additional property stating that the moments of order 0 and 1 of u and Π ∇ u coincide was added by B. Ahmad [44] to the space V 1 (Ω), resulting in the projection operator of Π 0 u being the same as the projection operator of Π ∇ u. For simplicity, Π P u is adopted to express Π ∇ u and Π 0 u uniformly. In this study, Π P u is deduced using an approach different from the standard procedure [46].
Beforehand, two notations are defined. For a function w, the average of the values it assumes at the vertices of Ω is denoted by w * : and the volume average of w over Ω is denoted by 〈w〉: where S is the area of Ω. P 1 (Ω) is the space of linear displacements over Ω, which can be split into two displacement spaces due to rigid motions and constant strain. We denote the space regarding rigid motions by R(Ω) and the space concerning constant strain by C(Ω): Obviously, P 1 (Ω) is the summation of R(Ω) and C(Ω). Both R(Ω) and C(Ω) are subspaces of V 1 (Ω). Two projection maps, Π R : V 1 (Ω)→R(Ω) and Π C : V 1 (Ω)→C(Ω), are defined to extract the rigid motions and constant strain of u∈V 1 (Ω). Considering that elements of C(Ω) should contain no rigid motion and that elements of R(Ω) should contain no constant strain, the following orthogonality conditions are imposed on the two maps: Gain et al. [48] proved that projection maps Π R and Π C satisfying the above properties can be given by and Due to the orthogonality conditions in Equations (28) and (29), we have It is noticed that and that ∇p in Equation (22) contains only constant elements. Obviously, the projection map Π P u given by Equation (32) satisfies Equation (22). Using the integration by parts, the areal integrals in Equation (32) can easily be converted to the boundary integrals, and then, it can be computed precisely because u = (u, v) is linear on each edge of Ω. Denoting the length of edge i connecting vertex i−1 Materials 2021, 14, 1252 9 of 24 and vertex i as l i , as Figure 4 shows, and the outward unit normal vector on edge i as n i = (n ix , n iy ) T , ∇u can be computed as Materials 2021, 14, x FOR PEER REVIEW 9 of 25 (u, v) is linear on each edge of Ω. Denoting the length of edge i connecting vertex i−1 and vertex i as li, as Figure 4 shows, and the outward unit normal vector on edge i as ni = (nix, niy) T , ⟨ u⟩ can be computed as l n l n l n l n u u u x x l n l n l n l n n S y y v v v For ease of subsequent performance, another matrix expression of ΠPu is given as in which and By substituting Equation (34) into Equation (32), we have For ease of subsequent performance, another matrix expression of Π P u is given as in which and . . . l n n ny +l 1 n 1y −4S l n n nx +l 1 n 1x 4S l 1 n 1x +l 2 n 2x 2S 0 l 2 n 2x +l 3 n 3x 2S 0 . . . l n n nx +l 1 n 1x 2S 0 0 l 1 n 1y +l 2 n 2y 2S 0 l 2 n 2y +l 3 n 3y 2S . . . 0 l n n ny +l 1 n 1y 2S l 1 n 1y +l 2 n 2y 2S l 1 n 1x +l 2 n 2x 2S l 2 n 2y +l 3 n 3y 2S l 2 n 2x +l 3 n 3x 2S . . . l n n ny +l 1 n 1y 2S l n n nx +l 1 n 1x 2S Due to the orthogonality conditions in Equations (28) and (29), the constant strain tensor ε(Π P u) can be computed as follows: and a matrix expression of the constant strain vector ε = (ε x , ε y , γ xy ) can be given as . . . 0 l n n ny +l 1 n 1y 2S l 1 n 1y +l 2 n 2y 2S l 1 n 1x +l 2 n 2x 2S l 2 n 2y +l 3 n 3y 2S l 2 n 2x +l 3 n 3x 2S . . . l n n ny +l 1 n 1y 2S l n n nx +l 1 n 1x 2S Obviously, P c is actually composed by the latter three rows in H. An element u∈V 1 (Ω) can be decomposed into Π P u and the residual item u−Π P u. Then, a bilinear form from V 1 (Ω) × V 1 (Ω) to R can be analyzed without the explicit expression of u, which will be demonstrated in the next section.

Computation of α(u, v) and f(v)
For simplicity of expression, α (u, v) in Equation (4) and f (v) in Equation (8) are rewritten as in which and It is noticed that, in this study, V 0 is considered an element in V 1 (Ω) to estimate V 0 within the block because V 0 at the block vertices are directly related to u. Consequently, α V is a bilinear form. V 0 can also be evaluated using the rigid body motion and the constant strain associated with Π P u. That will induce V 0 to become an element in P 1 (Ω) and α V to become a linear form.
The bilinear forms are formulated first. Due to the orthogonality conditions (22) and (23), any bilinear form α(u, v) regarding two elements u, v∈V 1 (Ω) can be rewritten as On the right side of Equation (45), the first term named the consistency term can be precisely calculated while the second term called the stabilization term represents the contribution of the residual item u−Π P u to the bilinear form. The stabilization term cannot produce an exact solution, so the second term in VEM is usually denoted as s(u−Π P u, v−Π P v) to indicate that it is an approximation. The standard scheme of s(u−Π P u, v−Π P v) was provided in classical VEM literature [46], but it is quite free in the construction of s(u−Π P u, v−Π P v) in practice [48,49]. The stabilization term proposed by Veiga [46] is adopted in this study: where η is a positive parameter ensuring the right scaling of the bilinear form assigned to the residual item. Using Equation (45), the bilinear form α E is rewritten as follows: Substituting the matrix expression (40) regarding ε(Π P u) into Equation (47), the consistency term is in which D is the elastic matrix The stabilization term in Equation (47) is where I is a 2n × 2n unit matrix and P u is a 2n × 2n constant matrix: In the view that s E is the approximation of α E , the positive parameters η E can be decided upon by requiring that s E and α E are comparable. The consistency term α E (Π P u, Π P u) can be estimated using s E as follows: The exact solution of α E (Π P u, Π P u) is given in Equation (48). Equating the traces of the two matrices, η E is determined: Combining Equations (48) and (49), we have Mimicking what we did for the bilinear form α E , the bilinear forms α M and α V are formulated as and respectively. Here, P m is a 2n × 2n constant matrix as and The integrations in Equation (56) can be exactly calculated by using the simplex integration. Now, we have successfully formulated all bilinear forms in Equation (43).
Then, the linear forms f σ , f b , and f t in Equation (44) are investigated. Obviously, f t can be computed directly. Because u is linear on each edge of Ω, a definite shape function matrix N(x) can be easily determined for a point x on ∂Ω to compute u(x) using Equation (6). Thus, we have Only the constant stress associated with Π P u can be valued in the first-order VEM, so σ 0 is a constant stress in this study. Using the integration by parts, the linear form f σ can be processed as Treating the constant body force b(x, y) as an element in P 1 (Ω), we have in which L b is a 2 × 2 n matrix as follows: and So far, all bilinear and linear forms in Equations (4) and (8) were formulated. The energy induced by the contact constraints are investigated in next section.

Computation of J'(v) Due to the Contact Constraints
In the block system, no tension and interpenetration is allowable when the collision of blocks occurs. Three contact states, "open", "lock", and "slide", are defined in DDA. By adopting stiff springs, different contact constraints are imposed according to the state.
A typical contact pair that is composed by a vertex i of Ω B and an entrance line jk of Ω A is plotted in Figure 5a If d n has a negative value, the contact pair must be "open" and no contact constraint applies, which results in no spring being added into the block system. Otherwise, a normal stiffness spring is added to oppose the penetration, which results in a deformation energy J n . Suppose that the two blocks contact each other along a coincident edge. The Coulomb friction law is adopted to determine whether the tangential relative motion is permitted. If the tangential relative motion is inadmissible, the contact pair is "lock" and a stiffness spring along the tangential direction is introduced to resist the tangential relative motion, leading to a deformation energy J τ . Otherwise, the contact pair is "slide" and a pair of friction forces is adopted along the tangential direction, inducing a potential energy J f . The allowable upper bound of the Coulomb friction law is considered the value of the friction forces. For the block system, J'(v) in Equation (1) can be obtained by collecting J n , J τ , and J f of all contact pairs. applies, which results in no spring being added into the block system. Otherwise, a normal stiffness spring is added to oppose the penetration, which results in a deformation energy J n . Suppose that the two blocks contact each other along a coincident edge. The Coulomb friction law is adopted to determine whether the tangential relative motion is permitted. If the tangential relative motion is inadmissible, the contact pair is "lock" and a stiffness spring along the tangential direction is introduced to resist the tangential relative motion, leading to a deformation energy J τ . Otherwise, the contact pair is "slide" and a pair of friction forces is adopted along the tangential direction, inducing a potential energy J f . The allowable upper bound of the Coulomb friction law is considered the value of the friction forces. For the block system, J'(v) in Equation (1) can be obtained by collecting J n , J τ , and J f of all contact pairs.
In practice, the states of all contact pairs have to be supposed in advance, and then, the "open-close" iteration is executed to guarantee the validity of the states. First, assemble the global equations based on the assumed contact states and obtain an interim displacement solution. Next, check the states of the contact pairs by using the interim displacement solution and alter the previously supposed states. Modify the contact constraints concerning every contact pair. Then, renew and resolve the global equations to obtain an updated interim displacement solution. Finally, stop the iteration when the actual states of all contact pairs agree with the supposed states and acquire the ultimate displacement solution.  In practice, the states of all contact pairs have to be supposed in advance, and then, the "open-close" iteration is executed to guarantee the validity of the states. First, assemble the global equations based on the assumed contact states and obtain an interim displacement solution. Next, check the states of the contact pairs by using the interim displacement solution and alter the previously supposed states. Modify the contact constraints concerning every contact pair. Then, renew and resolve the global equations to obtain an updated interim displacement solution. Finally, stop the iteration when the actual states of all contact pairs agree with the supposed states and acquire the ultimate displacement solution.
Firstly, considering that the displacement should be small in a time step, the normal penetration d n has an approximation (61) by neglecting the second-order infinitesimal value: in which Given that the stiffness of the normal spring is ρ n , the normal spring resisting the normal penetration induces the deformation energy as follows: Minimization of J n causes three second-order vectors to be added to the global load vector and nine second-order square matrices to be added into the global stiffness matrix.
Then, similar to d n , the tangential relative displacement d τ can be formulated as Provided that the stiffness of the tangential spring is ρ τ , the tangential spring results in the deformation energy: Minimization of J τ also leads to three second-order vectors being added to the global load vector and nine second-order square matrices being added into the global stiffness matrix.
Finally, if the contact pair is "slide", a pair pf friction forces is imposed along the tangential direction instead of the tangential spring. The unit direction vector from point j to point k is The allowable upper bound of the Coulomb friction law is denoted by f up . For d τ in Figure 5b, the friction force f up τ acts on Ω B at point i and the friction force-f up τ acts on Ω A at point o. The potential energy caused by the pair of friction force can be formulated as Minimization of J f causes three second-order vectors to be added to the global load vector. It is noticed that the Coulomb model of friction assumes that the sliding frictional force is proportional to the normal contact force. Provided that the joints between blocks have the mechanical properties of the friction angle ϕ and the inner cohesion c, f up is equivalent to |ρ n d n |tanϕ + c. If Equation (64) is adopted to estimate f up , the frictional energy term will contribute not only to the global force vector but also to the stiffness matrix, leading to non-symmetry of the stiffness matrix. To avoid this issue, f up is estimated by using d n obtained in the previous "open-close" iteration. In the view that the variation in d n is small when the "open-close" iteration is about to converge, such a simplification is acceptable.
In the original DDA, the incremental vertex displacements are estimated using Equation (11), while in the proposed method, the incremental vertex displacements are selected as the degrees of freedom. Therefore, the contributions of the contact restraints to the global stiffness matrix and load vector have simpler but more precise expressions in the proposed method than in the original DDA [36].
Minimize the potential energy-induced contact pairs in the block system one by one, and the contribution of J'(v) to the global equations are obtained and Equation (9) is fulfilled.
When Equation (9) is fulfilled completely, the global equations are resolved as the traditional DDA. For one time step, the vertex coordinates are updated by adding their incremental displacements to the previous coordinates once the open-close iteration converges. In addition, the constant stress associated with accumulation of the constant strain in Π P u over the foregoing time steps is considered the initial stress for the following time step.

Numerical Examples
Some numerical examples are analyzed to test the proposed method in this section. To evaluate the correctness, precision, and efficiency of the new approach, DDA with the post-adjustment strategy and vertex displacements-based DDA based on Wachspress interpolation [36] are adopted to solve these examples. For an objective comparison on the computational efficiency, all computations are conducted on the same laptop with 8 GB of RAM and Intel Core i7-6500U CPU. When conducting contact computation, ρ n was prescribed as 100E and ρ τ was prescribed as 40E for all methods. Here, E is Young's modulus of the block.

Rotating Triangular Block Problem
As drawn in Figure 6a, a rotating triangular block problem was designed by the authors to investigate the contact states in the equilibrium equation and in the resulting configurations. There was a triangular block on a rectangular ramp. The block and the ramp had the same material properties: ρ = 2.8 g/cm 3 , E = 20 GPa, and Poisson's ratio ν = 0.25. The gravity force of this model was ignored, and a concentrated load P = 10 kN acted at point 5 along the horizontal direction. All vertices of the foundation were fixed as the boundary conditions. 10 2.0631 × 10 −4 Although the rotational angle r0 is very tiny in each time step, the difference is sometimes considerable in the two normal penetration values. Moreover, the normal penetration value is negative in the resulting configurations at steps 5, 6, and 7, which violates the actual contact conditions completely. Then, the example was reanalyzed using the proposed method. For the contact pair point 7-Edge 14, the normal penetration values are also given in Table 1. The results indicated that the two normal penetration values are identical when the incremental displacements at the vertices are directly selected as the degrees of freedom. Firstly, DDA with the post-adjustment strategy was adopted to simulate the dynamic behavior of the model by using the max step displacement ratio δ = 0.1 and ∆ = 0.01 s. Figure 6b illustrates the configurations at the end of the 6, 9, 13, and 19 time steps. The triangular block rotated around point 7 for the eccentric moment from P. In the model, edge 14 had no rotation for its two end points that were both fixed. For the contact pair point 7-edge 14, we noted the normal penetration value in the equilibrium equations and measured the normal penetration value in the updated configurations. The results are listed in Table 1.
Although the rotational angle r 0 is very tiny in each time step, the difference is sometimes considerable in the two normal penetration values. Moreover, the normal penetration value is negative in the resulting configurations at steps 5, 6, and 7, which violates the actual contact conditions completely.
Then, the example was reanalyzed using the proposed method. For the contact pair point 7-Edge 14, the normal penetration values are also given in Table 1. The results indicated that the two normal penetration values are identical when the incremental displacements at the vertices are directly selected as the degrees of freedom.
Finally, vertex displacement-based DDA based on Wachspress interpolation [36] was used to solve this example and produced the same results as the proposed method. In this example, the time consumed by DDA with the post-adjustment strategy was 0.034 s, vertex displacement-based DDA in [36] spent 0.038 s, and the proposed method consumed 0.037 s.

Sliding Problem
There was a rectangular block on a triangular ramp, as Figure 7a shows. The top surface of the ramp had a slope angle of 45 • . The block and ramp were prescribed to have the same mechanical parameters: ρ = 2.5 × 10 3 kg/m 3 , E = 35 MPa, and ν = 0.3. The external load only came from the gravitational force with gravity acceleration g = 9.8 m/s 2 . No inner cohesion was involved in this model. All vertices of the ramp were fixed as the boundary conditions.  There was a rectangular block on a triangular ramp, as Figure 7a shows. The top surface of the ramp had a slope angle of 45°. The block and ramp were prescribed to have the same mechanical parameters: ρ = 2.5 × 10 3 kg/m 3 , E = 35 MPa, and ν = 0.3. The external load only came from the gravitational force with gravity acceleration g = 9.8 m/s 2 . No inner cohesion was involved in this model. All vertices of the ramp were fixed as the boundary conditions. The block slides along the ramp if the friction angle of the interface φ is less than 45°, and its displacement is  Figure 7b.
Above all, the relative errors of the two methods both appear to be correlated with the friction angle. Large friction angles always induce larger errors than small angles. The reason behind that is the approximation strategy of the contact force adopted in DDA. For a new contact pair, the proposed method initializes the contact force to 0 as the classical DDA and approximates the exact value according to the resulting interpenetration or inconsistent tangential motion between blocks.
Then, the relative errors of the two methods both rapidly drop off as the time steps increase. In the cases φ = 15° and φ = 30°, the displacement differences in the two methods are quite small. Both methods induced a relative error decreasing from 0.65% to 0.06% when φ = 15° and from 2.4% to 0.2% when φ = 30°, making it difficult to distinguish their relative error curves in Figure 7b. When φ = 40°, the original DDA results in a relative error of 55% at the first time step and a relative error of 2.9% at the end while the new method leads to a relative error of 37% at the first time step and a relative error of 2.0% at the end. Therefore, the proposed method provides higher accuracy results for this problem than the original DDA.  Figure 7b.
Above all, the relative errors of the two methods both appear to be correlated with the friction angle. Large friction angles always induce larger errors than small angles. The reason behind that is the approximation strategy of the contact force adopted in DDA. For a new contact pair, the proposed method initializes the contact force to 0 as the classical DDA and approximates the exact value according to the resulting interpenetration or inconsistent tangential motion between blocks.
Then, the relative errors of the two methods both rapidly drop off as the time steps increase. In the cases ϕ = 15 • and ϕ = 30 • , the displacement differences in the two methods are quite small. Both methods induced a relative error decreasing from 0.65% to 0.06% when ϕ = 15 • and from 2.4% to 0.2% when ϕ = 30 • , making it difficult to distinguish their relative error curves in Figure 7b. When ϕ = 40 • , the original DDA results in a relative error of 55% at the first time step and a relative error of 2.9% at the end while the new method leads to a relative error of 37% at the first time step and a relative error of 2.0% at the end. Therefore, the proposed method provides higher accuracy results for this problem than the original DDA.
As in the original DDA, the constant stress is estimated using the proposed method because only the constant strain associated with Π P u can be determined within the block. To investigate the difference in stress results of DDA with the post-adjustment strategy and the proposed method, the sliding problem was reanalyzed. Because the stress is unstable in a dynamic analysis, five friction angles, 50 • , 55 • , 60 • , 65 • , and 70 • were adopted and static analysis with 50 time steps was executed to obtain a stable stress result. Firstly, little difference occurs in the resulting stress values of the two methods. When ϕ = 50 • , the maximal principal stress σ 1 = 0.50 kPa and the minor principal stress σ 2 = −6.91 kPa in the original DDA and σ 1 = 0.63 kPa and σ 2 = −6.66 kPa in the proposed method. When ϕ = 70 • , σ 1 = 2.22 kPa and σ 2 = −5.97 kPa in the original DDA and σ 1 = 2.22 kPa and σ 2 = −5.91 kPa in the proposed method. Then, the minimum principal stress direction appears to increase gently with the friction angle increasing, as shown in Figure 8a, and it stabilizes when the friction angle increases to 65 • . The variations in the minimum principal stress direction with the friction angles are plotted for the two methods in Figure 8b. The proposed method has a gentler direction for the minimum principal stress. The stable direction for the minimum principal stress is −76.93 • in the original DDA and −76.60 • in the proposed method. methods in Figure 8b. The proposed method has a gentler direction for the minimum principal stress. The stable direction for the minimum principal stress is −76.93° in the original DDA and −76.60° in the proposed method.
For this problem, the calculation times consumed by the three methods appear to have a clear gap. Taking the case φ = 30° as an example, DDA with the post-adjustment strategy took 0.041 s. The time consumed by vertex displacement-based DDA in [36] is 0.108 s. This is because numerical integrations is necessary in computation for a rectangular block when using vertex displacement-based DDA based on Wachspress interpolation. However, only 0.042 s is consumed by vertex displacement-based DDA using VEM.

Surrounding Rock Problem
The stability assessment of a surrounding rock is an important problem for the design and construction of underground engineering. Figure 9 describes a surrounding rock model composed of 36 rock blocks. Both convex and concave blocks were included in this model. All blocks were assumed to have the same mechanical characteristics: ρ = 2.8 g/cm 3   For this problem, the calculation times consumed by the three methods appear to have a clear gap. Taking the case ϕ = 30 • as an example, DDA with the post-adjustment strategy took 0.041 s. The time consumed by vertex displacement-based DDA in [36] is 0.108 s. This is because numerical integrations is necessary in computation for a rectangular block when using vertex displacement-based DDA based on Wachspress interpolation. However, only 0.042 s is consumed by vertex displacement-based DDA using VEM.

Surrounding Rock Problem
The stability assessment of a surrounding rock is an important problem for the design and construction of underground engineering. Figure 9 describes a surrounding rock model composed of 36 rock blocks. Both convex and concave blocks were included in this model. All blocks were assumed to have the same mechanical characteristics: ρ = 2.8 g/cm 3  Using δ = 0.005 and a self-adjusting Δ from the code, the proposed method was employed to conduct a static analysis of the problem. Taking the friction angle φ = 30° for the joints between blocks, Figure 10a illustrates the ultimate configurations and main stresses after 100 time steps. The proposed method arrived at the conclusion that the model is stable. The stress level decreases from the outer layer to the inner layer. At the top and the bottom, the inner layer suffers the largest horizontal compress effect. The top two rocks in the inner layer suffer σ1 = −0.4 kPa and σ2 = −109.7 kPa, which induces a minimum principal stress direction of 6.7°. The bottom two rocks in the inner layer suffer Using δ = 0.005 and a self-adjusting ∆ from the code, the proposed method was employed to conduct a static analysis of the problem. Taking the friction angle ϕ = 30 • for the joints between blocks, Figure 10a illustrates the ultimate configurations and main stresses after 100 time steps. The proposed method arrived at the conclusion that the model is stable. The stress level decreases from the outer layer to the inner layer. At the top and the bottom, the inner layer suffers the largest horizontal compress effect. The top two rocks in the inner layer suffer σ 1 = −0.4 kPa and σ 2 = −109.7 kPa, which induces a minimum principal stress direction of 6.7 • . The bottom two rocks in the inner layer suffer σ 1 = −15.6 kPa and σ 2 = −138.6 kPa, and the direction of the minimum principal stress is 11.4 • . Figure 9. A surrounding rock model.
Using δ = 0.005 and a self-adjusting Δ from the code, the proposed method was employed to conduct a static analysis of the problem. Taking the friction angle φ = 30° for the joints between blocks, Figure 10a illustrates the ultimate configurations and main stresses after 100 time steps. The proposed method arrived at the conclusion that the model is stable. The stress level decreases from the outer layer to the inner layer. At the top and the bottom, the inner layer suffers the largest horizontal compress effect. The top two rocks in the inner layer suffer σ1 = −0.4 kPa and σ2 = −109.7 kPa, which induces a minimum principal stress direction of 6.7°. The bottom two rocks in the inner layer suffer σ1 = −15.6 kPa and σ2 = −138.6 kPa, and the direction of the minimum principal stress is 11.4°. Considering the mechanical property of the joints may be damaged by the blasting excavation effect, friction was neglected and this model was reanalyzed. Figure 10b plots the results after 100 time steps, which adopts the same tags for stress vectors as Figure  10a. Because of the arch effect, the surrounding rock is also stable and the stress level still decreases from the outer to the inner layers. In contrast, the stress values appear to increase in the case of no friction. The top two rocks in the inner layer suffer σ1 = −6.2 kPa and σ2 = −117.6 kPa, and the minimum principal stress direction is 9.2°. The bottom two rocks in the inner layer suffer σ1 = −24.8 kPa and σ2 = −227.5 kPa, which has a minimum principal stress direction of 8.4°. The reason for this is that the stability of the arch only relies on the compression between blocks when the joints friction effect is absent. Considering the mechanical property of the joints may be damaged by the blasting excavation effect, friction was neglected and this model was reanalyzed. Figure 10b plots the results after 100 time steps, which adopts the same tags for stress vectors as Figure 10a. Because of the arch effect, the surrounding rock is also stable and the stress level still decreases from the outer to the inner layers. In contrast, the stress values appear to increase in the case of no friction. The top two rocks in the inner layer suffer σ 1 = −6.2 kPa and σ 2 = −117.6 kPa, and the minimum principal stress direction is 9.2 • . The bottom two rocks in the inner layer suffer σ 1 = −24.8 kPa and σ 2 = −227.5 kPa, which has a minimum principal stress direction of 8.4 • . The reason for this is that the stability of the arch only relies on the compression between blocks when the joints friction effect is absent.
The stability of the surrounding rock in this problem was also verified by using DDA with the post-adjustment strategy. Vertex displacement-based DDA based on Wachspress interpolation [36] fails to solve this problem because Wachspress interpolation does not apply for concave blocks. For this problem, more calculation time is consumed by the proposed method than DDA with the post-adjustment strategy, which can be attributed to the fact that the proposed method almost doubles DDA with the post-adjustment strategy in the dimension of the total degrees of freedom. Taking the case with no friction as an example, DDA with the post-adjustment strategy took 0.891 s and the proposed method took 1.117 s.

Block Wall Failure Problem
As Figure 11 shows, a wall containing 70 blocks was analyzed to test the capacity of the proposed method. The wall had a width of 270 m and a height of 180 m. The mechanical parameters of the blocks were given as ρ = 2.0 × 10 3 kg/m 3 , E = 500 MPa, and ν = 0.25. The external load only included the gravitational force with g = 10 m/s 2 . No friction occurred between the blocks. The outer vertices of two base blocks were fixed as the boundary conditions.
Taking ∆ = 0.012 s and δ = 0.01, the dynamic behavior of the model was simulated. Figure 12 demonstrates the resulted configurations and main stresses after 7 s.

Block Wall Failure Problem
As Figure 11 shows, a wall containing 70 blocks was analyzed to test the capacity of the proposed method. The wall had a width of 270 m and a height of 180 m. The mechanical parameters of the blocks were given as ρ = 2.0 × 10 3 kg/m 3 , E = 500 MPa, and ν = 0.25. The external load only included the gravitational force with g = 10 m/s 2 . No friction occurred between the blocks. The outer vertices of two base blocks were fixed as the boundary conditions. Figure 11. The block wall model.
Taking Δ = 0.012 s and δ = 0.01, the dynamic behavior of the model was simulated. Figure 12 demonstrates the resulted configurations and main stresses after 7 s. First, at the bottom, the central blocks descend because of the absence of supports. Next, the middle blocks subside one layer after another. Simultaneously, the movement of the middle blocks leads to rotation and inclination of the left and right blocks. Submitting to the gravitational force of the upper blocks, the bottom two blocks become deformed. After 7 s elapsed, the top boundary of the model appears to be a gentle curve. The stress value of the blocks in the middle is rather small compared to the blocks in other regions, and the two bottom blocks suffer the largest stress. In the bottom right block, the main stresses are σ1 = 560 kPa and σ2 = −880 kPa and the direction of the minimum principal stress is −48°.
With Δ and δ unchanged, the dynamic behavior of the example was reanalyzed using DDA with the post-adjustment strategy. The configurations and main stress vectors after 7 s elapsed are shown in Figure 13. First, at the bottom, the central blocks descend because of the absence of supports. Next, the middle blocks subside one layer after another. Simultaneously, the movement of the middle blocks leads to rotation and inclination of the left and right blocks. Submitting to the gravitational force of the upper blocks, the bottom two blocks become deformed. After 7 s elapsed, the top boundary of the model appears to be a gentle curve. The stress value of the blocks in the middle is rather small compared to the blocks in other regions, and the two bottom blocks suffer the largest stress. In the bottom right block, the main stresses are σ 1 = 560 kPa and σ 2 = −880 kPa and the direction of the minimum principal stress is −48 • .
With ∆ and δ unchanged, the dynamic behavior of the example was reanalyzed using DDA with the post-adjustment strategy. The configurations and main stress vectors after 7 s elapsed are shown in Figure 13.
The stress value of the blocks in the middle is rather small compared to the blocks in other regions, and the two bottom blocks suffer the largest stress. In the bottom right block, the main stresses are σ1 = 560 kPa and σ2 = −880 kPa and the direction of the minimum principal stress is −48°.
With Δ and δ unchanged, the dynamic behavior of the example was reanalyzed using DDA with the post-adjustment strategy. The configurations and main stress vectors after 7 s elapsed are shown in Figure 13.  The difference in the results of the proposed method and DDA with the post-adjustment strategy is distinct. The two bottom blocks suffer larger stresses in DDA with the postadjustment strategy than in the proposed method. In the bottom right block, σ 1 = 2115 kPa and σ 2 = −2205 kPa, and the minimum principal stress has a direction of −44 • . Rotation and inclination of the blocks in the left and right regions are smaller compared to the results of the proposed method. The difference in the displacement mode of the two methods has great effects on the difference between the results. The linear displacement function was adopted directly in DDA with the post-adjustment strategy. However, in the case when the vertex number was n > 3, the contribution of the higher-order displacement was added to the proposed method. Therefore, the proposed approach provides a block of n > 3 vertices with larger deformability than DDA with the post-adjustment strategy. Once the deformation of this block is not restrained by the adjacent blocks, a larger deformation occurs in the proposed method.
For this problem, DDA with the post-adjustment strategy took 1.606 s, the new method took 1.943 s, and vertex displacement-based DDA based on Wachspress interpolation [36] took 4.762 s.
The results of Examples 1 and 2 verified that the proposed method has a higher accuracy in the contact computation than DDA with the post-adjustment strategy. Previous studies to construct vertex displacement-based DDA focused on governing the displacement of a block using the displacement functions in PFEM, e.g., Wachspress interpolation [36]. The issue that Wachspress interpolation does not apply for concave blocks can be remedied by introducing those generalized barycentric coordinates well-defined in the concave domain. If the displacement functions in PFEM are adopted in the construction of vertex displacement-based DDA, the numerical integrations have to be executed for a block with more than three vertices, inducing lower computational efficiency for large-scale problems. Relatively speaking, the proposed method may be a better choice when taking the displacements at the vertices as the degrees of freedom. In the view that the degrees of freedom have a larger dimension in the proposed method than in the original, the new method has a lower computational efficiency than the original DDA. However, a 2% to 20% increase in time consumption is achieved, which is acceptable.

Conclusions
In this study, the increments of displacements at block vertices in one time step were selected as the new degrees of freedom for a block. The virtual element spaces were defined for a block and the block system. In the framework of VEM theory, the total potential energy of the DDA block system was estimated and minimized to obtain the global equilibrium equation and load vector. The proposed method can be considered an improvement of DDA by reformulating the block displacement using VEM and can be regarded as an attempt to model the dynamic behavior of the block system using VEM coupled with the contact theory in DDA.
The proposed method has some advantages. On the one hand, the approximation of cos r 0 and sin r 0 in the displacement function of the original DDA is avoided in the proposed method. The degrees of freedom are directly added to the former coordinates in the renewal of the vertex coordinates, which induces a simpler expression regarding the contact restraint impositions. On the other hand, the contact detection and the contact theory can still be adopted as in the original DDA. Meanwhile, most of the numerical integrations are executed along the block boundary when calculating the global stiffness matrix and load vector. Only numerical integrations within the block can be conducted using the simplex integration as in the classical DDA.
Additionally, V 1 (Ω) is taken as the virtual element space for a block in this study. If the stress variability within the block is expected, a higher-order virtual element space instead of V 1 (Ω), e.g., V 2 (Ω) and V 3 (Ω), can be adopted for a block. When using a higher-order virtual element space, the degrees of freedom should be augmented; refer to [46]. Although this study is limited to 2D analysis, the method to improve DDA by reformulating the block displacement using VEM can be mimicked in a 3D case.