Next Article in Journal
Fatigue Life of 7005 Aluminum Alloy Cruciform Joint Considering Welding Residual Stress
Next Article in Special Issue
Concurrent Lamination and Tapering Optimization of Cantilever Composite Plates under Shear
Previous Article in Journal
Short-Term Impact of AC Harmonics on Aging of NiMH Batteries for Grid Storage Applications
Previous Article in Special Issue
Homogenization of Radial Temperature by a Tungsten Sink in Sublimation Growth of 45 mm AlN Single Crystal
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Vertex Displacement-Based Discontinuous Deformation Analysis Using Virtual Element Method

1
State Key Laboratory of Geomechanics and Geotechnical Engineering, Institute of Rock and Soil Mechanics, Chinese Academy of Sciences, Wuhan 430071, China
2
State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin, China Institute of Water Resources and Hydropower Research, Beijing 100038, China
3
Key Laboratory of Geological Hazards on Three Gorges Reservoir Area, Ministry of Education, China Three Gorges University, Yichang 443002, China
*
Author to whom correspondence should be addressed.
Submission received: 31 December 2020 / Revised: 17 February 2021 / Accepted: 2 March 2021 / Published: 6 March 2021
(This article belongs to the Special Issue Numerical Methods and Optimization of Structures: FEM)

Abstract

:
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.

1. 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 r0 = r0 and cos r0 = 1 for a small rotation angle r0, 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 V1(Ω) is adopted to describe the displacements of points in the block, and the projector ΠPu from V1(Ω) on the linear displacement space P1(Ω) 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 ΠPu and an approximation of uPu. 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 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.

2. 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 Γ u I and the traction boundary Γ t I . The displacement on Γ u I is prescribed as ûI, and the surface traction on Γ t I is denoted by tI. 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.
Within a time step, DDA requires the incremental displacements u = (u1, u2, …, unb) to be the unique minimizer of the total potential energy of the block system:
u = arg min v V J ( v ) ,   J ( v ) = 1 2 α ( v , v ) f ( v ) + J ( v )
in which V is the incremental displacement space of the block system:
  V = V 1 × V 2 × × V n b ,   with   V I = { v I [ H 1 ( Ω I ) ] 2 : v I | Γ u I = u ^ I }
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 uIVI and uLVL.
The bilinear form α(u, v) and the linear form f(v) are computed for u∈V and v∈V by
α ( u , v ) = I = 1 n b α ( u I , v I ) , f ( v ) = I = 1 n b f ( v I )
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 uV and vV represents the energy due to elastic deformation and the inertial force:
α ( u , v ) = Ω σ ( u ) : ε ( v )   d x + Ω ρ u ¨ v   d x
where ε is the incremental constant strain decided by u:
ε ( u ) = 1 2 ( u + T u )
in which ∇ is the gradient operator. σ is the incremental Cauchy stress related to ε by the constitutive equation:
σ ( u ) = D [ ε ( u ) ]
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 V0 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
u ¨ = 2 Δ 2 u 2 Δ V 0
The linear form f(v) is defined by
f ( v ) = Ω b v d x Ω σ 0 : ε ( v ) d x + Ω t v d s
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.
[ K 11 K 12 K 13 K 1 ( n b ) K 21 K 22 K 23 K 2 ( n b ) K 31 K 32 K 33 K 3 ( n b ) K ( n b ) 1 K ( n b ) 2 K ( n b ) 3 K ( n b ) ( n b ) ] [ d 1 d 2 d 3 d n b ] = [ F 1 F 2 F 3 F n b ]
in which dI is the degrees of freedom concerning block I. FI is the generalized load vector with the same dimension of dI. Denoting the dimension of dI with dim(dI), KII is a dim(dI) × dim(dI) matrix and KIL(I≠L) is a dim(dI) × dim(dL) 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, KIL is zero.

3. Demerits Caused by the Original Degrees of Freedom and Previous Attempts to Construct a Vertex Displacement-Based DDA

3.1. Demerits Caused by the Original Degrees of Freedom in DDA

For block Ω, the traditional DDA defines the degrees of freedom d within one time step:
d = ( u 0 v 0 r 0 ε x ε y γ x y ) T
where u0 and v0 represent the increments of the rigid horizontal and vertical translations, respectively; r0 denotes the increment of the rigid rotation angle around the central point (x0, y0) 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 = T d  
and T is the translation matrix as follows:
T = [ 1 0 ( y y o ) ( x x o ) 0 0.5 ( y y o ) 0 1 ( x x o ) 0 ( y y o ) 0.5 ( x x o ) ]
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
u = u 0 + ( x x 0 ) ( cos r 0 1 ) ( y y 0 ) sin r 0 + ( x x 0 ) ε x + 0.5 ( y y 0 ) γ x y v = v 0 + ( x x 0 ) sin r 0 + ( y y 0 ) ( cos r 0 1 ) + ( y y 0 ) ε y + 0.5 ( x x 0 ) γ x y
Comparing Equation (13) with Equation (11), the approximations of cos r0 = 1 and sin r0 = r0 are adopted for the small rotation angle r0 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.

3.2. 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 ui= (ui, vi)T at the vertex xi constitute the new degrees of freedom:
d ¯ = ( u 1 v 1 u 2 v 2 u n v n ) T
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
u = ( u v ) = ( N 1 0 N 2 0 N n 0 0 N 1 0 N 2 0 N n ) d ¯ = N d ¯
where Ni is the shape function in terms of the vertex xi. Because the block in DDA might have more than four vertices, the shape functions developed in PFEM were introduced to construct Ni [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 Ni to develop an effective PFEM [38]. As shown in Figure 2, a block Ω has n vertices located at xi and its boundary is composed of n straight edges. The red edges with vertex xi are denoted by Γi. Ni should be satisfied.
(1) Partition of unity, boundedness, and nonnegativity:
i = 1 n N i ( x ) = 1 , 0 N i ( x ) 1
(2) Kronecker-delta property:
  N i ( x j ) = δ i j
(3) Linear completeness:
i = 1 n N i ( x )   x i = x
(4) With C being within block Ω while C0 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.
N i ( x ) = A ( x i 1 , x i , x i + 1 ) ( Π k i , i 1 A ( x , x k , x k + 1 ) ) j = 1 n A ( x j 1 , x j , x j + 1 ) ( Π k j , j 1 A ( x , x k , x k + 1 ) )
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 the interpolation function Ni in PFEM. When Equation (15) is adopted to govern the block displacement, some advantages of the original DDA can compromised.

4. Vertex Displacement-Based DDA Using VEM

4.1. 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:
  V k ( Ω ) = { v [ H 1 ( Ω ) C 0 ( Ω ) ] 2 : Δ v [ P k 2 ( Ω ) ] 2 , v | Γ [ P k ( Γ ) ] 2 Γ Ω } ,
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(Ω)]2Vk(Ω). For simplicity, [Pk(Ω)]2 is denoted as Pk(Ω).
We take k = 1 in this study. There are three main reasons for this choice: (i) The degrees of freedom for V1(Ω) 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 V1(Ω). (iii) In the view that the original DDA adopted a complete first-order displacement function, the degree of accuracy is not threatened by using V1(Ω). Therefore, an individual virtual element space V1(Ω) 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 V1(Ω) is 2n, where n denotes the number of vertices for the block.
Provided that there are nb blocks in the block system, the total virtual element space can be given as
  V 1 = V 1 1 ( Ω 1 ) × V 1 2 ( Ω 2 ) × × V 1 n b ( Ω n b ) ,
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 V1.

4.2. The Projection Operator ΠPu: V1(Ω)→P1(Ω)

To minimize the total potential energy of the block system, the computation of α (u, v) from V1(Ω) × V1(Ω) to R defined in Equation (4) is inevitable, which can be accomplished in the framework of VEM theory. Referring to [47], two projectors Πu: V1(Ω)→P1(Ω) and Π0u: V1(Ω)→P1(Ω) are respectively defined using the orthogonality conditions:
Ω ( Π u ) p   dx = Ω u p   dx   p P 1 ( Ω ) ,
and
Ω ( Π 0 u ) p   dx = Ω u p   dx   p P 1 ( Ω ) .
where Π0u 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 V1(Ω), resulting in the projection operator of Π0u being the same as the projection operator of Πu. For simplicity, ΠPu is adopted to express Πu and Π0u uniformly. In this study, ΠPu 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:
w = 1 n i = 1 n w ( x i ) ,
and the volume average of w over Ω is denoted by ⟨w⟩:
w = 1 S Ω w   d x .
where S is the area of Ω.
P1(Ω) 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(Ω):
R ( Ω ) = { a + B R ( x x * ) : a R 2 , B R R 2 × 2 , B R T = B R } ,
C ( Ω ) = { B C ( x x * ) : B C R 2 × 2 , B C T = B C } .
Obviously, P1(Ω) is the summation of R(Ω) and C(Ω). Both R(Ω) and C(Ω) are subspaces of V1(Ω). Two projection maps, ΠR: V1(Ω)→R(Ω) and ΠC: V1(Ω)→C(Ω), are defined to extract the rigid motions and constant strain of uV1(Ω). 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:
Π R c = 0 ,   c C ( Ω ) ,
Π C r = 0 ,   r R ( Ω ) .
Gain et al. [48] proved that projection maps ΠR and ΠC satisfying the above properties can be given by
Π R u = u + ω ( u ) ( x x )   with   ω ( u ) = 1 2 ( u T u )
and
Π C u = ε ( u ) ( x x )  
Due to the orthogonality conditions in Equations (28) and (29), we have
Π P u = Π R u + Π C u = u + u ( x x ) .
It is noticed that
Ω ( Π P u )   d x = Ω u   d x =   Ω u   d x ,
and that ∇p in Equation (22) contains only constant elements. Obviously, the projection map ΠPu 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 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
u = 1 S Ω u   d x = 1 S [ Ω u n x   d s Ω u n y   d s Ω v n x   d s Ω v n y   d s ] = 1 S [ i = 1 n l i n i x + l i + 1 n ( i + 1 ) x 2 u i i = 1 n l i n i y + l i + 1 n ( i + 1 ) y 2 u i i = 1 n l i n i x + l i + 1 n ( i + 1 ) x 2 v i i = 1 n l i n i y + l i + 1 n ( i + 1 ) y 2 v i ]
By substituting Equation (34) into Equation (32), we have
Π P u = 1 n ( i = 1 n u i i = 1 n v i ) + 1 S [ i = 1 n l i n i x + l i + 1 n ( i + 1 ) x 2 u i i = 1 n l i n i y + l i + 1 n ( i + 1 ) y 2 u i i = 1 n l i n i x + l i + 1 n ( i + 1 ) x 2 v i i = 1 n l i n i y + l i + 1 n ( i + 1 ) y 2 v i ] ( x x y y )
For ease of subsequent performance, another matrix expression of ΠPu is given as
Π P u = T * H d ¯ ,
in which
T * = [ 1 0 ( y y * ) ( x x * ) 0 0.5 ( y y * ) 0 1 ( x x * ) 0 ( y y * ) 0.5 ( x x * ) ]
and
  H = ( 1 n 0 1 n 0 1 n 0 0 1 n 0 1 n 0 1 n l 1 n 1 y + l 2 n 2 y 4 S l 1 n 1 x + l 2 n 2 x 4 S l 2 n 2 y + l 3 n 3 y 4 S l 2 n 2 x + l 3 n 3 x 4 S l n n n y + l 1 n 1 y 4 S l n n n x + l 1 n 1 x 4 S l 1 n 1 x + l 2 n 2 x 2 S 0 l 2 n 2 x + l 3 n 3 x 2 S 0 l n n n x + l 1 n 1 x 2 S 0 0 l 1 n 1 y + l 2 n 2 y 2 S 0 l 2 n 2 y + l 3 n 3 y 2 S 0 l n n n y + l 1 n 1 y 2 S l 1 n 1 y + l 2 n 2 y 2 S l 1 n 1 x + l 2 n 2 x 2 S l 2 n 2 y + l 3 n 3 y 2 S l 2 n 2 x + l 3 n 3 x 2 S l n n n y + l 1 n 1 y 2 S l n n n x + l 1 n 1 x 2 S ) .
Due to the orthogonality conditions in Equations (28) and (29), the constant strain tensor εPu) can be computed as follows:
ε ( Π P u ) = ε ( Π C u ) = ε ( u ) = 1 S [ i = 1 n l i n i x + l i + 1 n ( i + 1 ) x 2 u i i = 1 n ( l i n i y + l i + 1 n ( i + 1 ) y 4 u i + + l i n i x + l i + 1 n ( i + 1 ) x 4 v i ) i = 1 n ( l i n i y + l i + 1 n ( i + 1 ) y 4 u i + + l i n i x + l i + 1 n ( i + 1 ) x 4 v i ) i = 1 n l i n i y + l i + 1 n ( i + 1 ) y 2 v i ]
and a matrix expression of the constant strain vector ε = (εx, εy, γxy) can be given as
ε ( Π P u ) = 1 S [ i = 1 n l i n i x + l i + 1 n ( i + 1 ) x 2 u i i = 1 n l i n i y + l i + 1 n ( i + 1 ) y 2 v i i = 1 n ( l i n i y + l i + 1 n ( i + 1 ) y 2 u i + + l i n i x + l i + 1 n ( i + 1 ) x 2 v i ) ] = P C d ¯
in which
P c = ( l 1 n 1 x + l 2 n 2 x 2 S 0 l 2 n 2 x + l 3 n 3 x 2 S 0 l n n n x + l 1 n 1 x 2 S 0 0 l 1 n 1 y + l 2 n 2 y 2 S 0 l 2 n 2 y + l 3 n 3 y 2 S 0 l n n n y + l 1 n 1 y 2 S l 1 n 1 y + l 2 n 2 y 2 S l 1 n 1 x + l 2 n 2 x 2 S l 2 n 2 y + l 3 n 3 y 2 S l 2 n 2 x + l 3 n 3 x 2 S l n n n y + l 1 n 1 y 2 S l n n n x + l 1 n 1 x 2 S ) .
Obviously, Pc is actually composed by the latter three rows in H.
An element uV1(Ω) can be decomposed into ΠPu and the residual item u−ΠPu. Then, a bilinear form from V1(Ω) × V1(Ω) to R can be analyzed without the explicit expression of u, which will be demonstrated in the next section.

4.3. 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
α ( u , v ) = α E + α M α V , f ( v ) = f b + f t f σ
in which
α E = Ω ε T ( v ) σ ( u )   d x ,   α M = 2 ρ Δ 2 Ω v T u   d x ,   α V = 2 ρ Δ Ω v T V 0   d x
and
f b = Ω v T b   d x ,   f t = Ω v T t d s ,   f σ = Ω ε T ( v ) σ 0 d x
It is noticed that, in this study, V0 is considered an element in V1(Ω) to estimate V0 within the block because V0 at the block vertices are directly related to u. Consequently, αV is a bilinear form. V0 can also be evaluated using the rigid body motion and the constant strain associated with ΠPu. That will induce V0 to become an element in P1(Ω) 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, vV1(Ω) can be rewritten as
α ( u , v ) = a ( Π P u , Π P v ) + a ( u Π P u , v Π P v ) .
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−ΠPu to the bilinear form. The stabilization term cannot produce an exact solution, so the second term in VEM is usually denoted as s(u−ΠPu, v−ΠPv) to indicate that it is an approximation. The standard scheme of s(u−ΠPu, v−ΠPv) was provided in classical VEM literature [46], but it is quite free in the construction of s(u−ΠPu, v−ΠPv) in practice [48,49]. The stabilization term proposed by Veiga [46] is adopted in this study:
s ( u Π P u , v Π P v ) = i = 1 n η [ u ( x i ) Π P u ( x i ) ] [ v ( x i ) Π P v ( x i ) ]
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:
α E ( u , u ) = α E ( Π P u , Π P u ) + s E ( u Π P u , u Π P u )
Substituting the matrix expression (40) regarding εPu) into Equation (47), the consistency term is
α E ( Π P u , Π P u ) = Ω ε T ( Π P u ) σ ( Π u )   d x = d ¯ T [ S P c T D P c ] d ¯
in which D is the elastic matrix
The stabilization term in Equation (47) is
s E ( u Π P u , u Π P u ) = i = 1 n η E [ u ( x i ) Π P u ( x i ) ] 2 = d ¯ T [ η E ( I P u ) T ( I P u ) ] d ¯
where I is a 2n × 2n unit matrix and Pu is a 2n × 2n constant matrix:
P u = [ T * ( x 1 ) T * ( x 2 ) T * ( x n ) ] H
In the view that sE is the approximation of αE, the positive parameters ηE can be decided upon by requiring that sE and αE are comparable. The consistency term αEPu, ΠPu) can be estimated using sE as follows:
s E ( Π P u , Π P u ) = d ¯ T [ η E P u T P u ] d ¯
The exact solution of αEPu, ΠPu) is given in Equation (48). Equating the traces of the two matrices, ηE is determined:
η E = t r a c e ( S P c T D P c ) t r a c e ( P u T P u )
Combining Equations (48) and (49), we have
α E ( u , u ) = d ¯ T [ S P c T D P c ] d ¯ + d ¯ T [ η E ( I P u ) T ( I P u ) ] d ¯
Mimicking what we did for the bilinear form αE, the bilinear forms αM and αV are formulated as
α M ( u , u ) = d ¯ T [ 2 ρ Δ 2 H T P m H ] d ¯ + d ¯ T [ η M ( I P u ) T ( I P u ) ] d ¯
and
α V ( u , V 0 ) = d ¯ T [ 2 ρ Δ H T P m H ] V 0 + d ¯ T [ η V ( I P u ) T ( I P u ) ] V 0 ,
respectively. Here, Pm is a 2n × 2n constant matrix as
P m = Ω ( T * ) T T * d x d y ,
η M = 2 ρ Δ 2 t r a c e ( H T P m H ) t r a c e ( P u T P u )
and
η V = 2 ρ Δ t r a c e ( H T P m H ) t r a c e ( P u T P u ) .
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σ, fb, and ft in Equation (44) are investigated. Obviously, ft 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
f t = d ¯ T Ω N T ( x ) t   d Γ t
Only the constant stress associated with ΠPu 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
f σ = Ω ε T ( u ) σ 0   d x = d ¯ T ( S P c T σ 0 )
Treating the constant body force b(x, y) as an element in P1(Ω), we have
f b = Ω u T b   d x = Ω ( Π P u ) T b   d x = d ¯ T [ H T Ω ( T * ) T   d x ] b = d ¯ T ( L b ) T b
in which Lb is a 2 × 2 n matrix as follows:
L b = [ L 1 b 0 L 2 b 0 L n b 0 0 L 1 b 0 L 2 b 0 L n b ]
and
L i b = S n + 0.5 ( l i n i x + l i + 1 n ( i + 1 ) x l i n i y + l i + 1 n ( i + 1 ) y ) ( x 0 x * y 0 y * )
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.

4.4. 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. Point o is the closest point from line jk to point i. At the start of a time step, the coordinates of points i, j, k, and o are (xi, yi), (xj, yj), (xk, yk), and (xo, yo), respectively. Within this time step, their displacements are (ui, vi), (uj, vj), (uk, vk), and (uo, vo), respectively. Then, the normal penetration dn and tangential relative displacement dτ in Figure 5b can be measured using their new coordinates when the time step terminates.
If dn 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 Jn. 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 Jf. 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 Jn, Jτ, and Jf 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.
Firstly, considering that the displacement should be small in a time step, the normal penetration dn has an approximation (61) by neglecting the second-order infinitesimal value:
d n = [ S 0 l + 1 l ( y j y k x k x j ) ( u i v i ) + 1 l ( y k y i x i x k ) ( u j v j ) + 1 l ( y i y j x j x i ) ( u k v k ) ]
in which
S 0 = | 1 x i y i 1 x j y j 1 x k y k | , l = ( x j x k ) 2 + ( y j y k ) 2
Given that the stiffness of the normal spring is ρn, the normal spring resisting the normal penetration induces the deformation energy as follows:
J n = 0.5 ρ n d n 2
Minimization of Jn 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 dn, the tangential relative displacement dτ can be formulated as
d τ = S ¯ 0 l + 1 l ( x k x j y k y j ) ( u i v i ) + μ l ( x j x k y j y k ) ( u j v j ) + 1 μ l ( x j x k y j y k ) ( u k v k )
where
S ¯ 0 = ( x k x j y k y j ) ( x i x 0 y i y 0 ) , μ = 1 l ( ( x 0 x k ) 2 + ( y 0 y k ) 2 )
Provided that the stiffness of the tangential spring is ρτ, the tangential spring results in the deformation energy:
J τ = 0.5 ρ τ d τ 2
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
τ = 1 l ( ( x k x j ) ( y k y j ) ) T
The allowable upper bound of the Coulomb friction law is denoted by fup. For dτ in Figure 5b, the friction force fupτ acts on ΩB at point i and the friction force—fupτ acts on ΩA at point o. The potential energy caused by the pair of friction force can be formulated as
J f = μ ( u j v j ) f u p τ ( 1 μ ) ( u k v k ) f u p τ + ( u i v i ) f u p τ
Minimization of Jf 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, fup is equivalent to |ρndn|tanφ + c. If Equation (64) is adopted to estimate fup, 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, fup is estimated by using dn obtained in the previous “open–close” iteration. In the view that the variation in dn 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 ΠPu over the foregoing time steps is considered the initial stress for the following time step.

5. 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.

5.1. 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/cm3, 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.
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 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.
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.

5.2. 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 × 103 kg/m3, E = 35 MPa, and ν = 0.3. The external load only came from the gravitational force with gravity acceleration g = 9.8 m/s2. 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
s = 0.5 g ( sin 45 o tan φ cos 45 o ) t 2
in which s denotes the sliding distance (m) of the rectangle and t is the elapsed time (s).
To appraise the accuracy of the new method, three values, 15°, 30°, and 40°, were assigned to the friction angle. Taking Δ = 0.02s and δ = 0.10, 35 time steps were computed. Referring to the analytical solutions of Equation (72), the relative errors of the proposed approach and the original DDA results were measured and are drawn in 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 ΠPu 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.
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.

5.3. 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/cm3, ν = 0.20, and E = 200 MPa. Only the gravitational force with g = 10 m/s2 acts on the model. The hoop joints in the surrounding rock divided the model into three layers. The outermost boundaries of this model were fixed as the displacement boundary conditions.
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.
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.

5.4. 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 × 103 kg/m3, E = 500 MPa, and ν = 0.25. The external load only included the gravitational force with g = 10 m/s2. 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.
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 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 post-adjustment 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.

6. 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 r0 and sin r0 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, V1(Ω) 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 V1(Ω), e.g., V2(Ω) and V3(Ω), 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.

Author Contributions

Conceptualization, methodology, H.L. and W.J.; software, G.S.; validation, formal analysis, investigation, resources, data curation, H.L., L.L., and W.J.; writing—original draft preparation, H.L.; writing—review and editing, W.J.; visualization, G.S.; funding acquisition, H.L. and W.J. All authors have read and agreed to the published version of the manuscript.

Funding

This research was funded by the Open Research Fund of State Key Laboratory of Simulation and Regulation of Water Cycle in River Basin (China Institute of Water Resources and Hydropower Research), grant number IWHR-SKL-202020, and by National Natural Science Funds, grant number 52079070 and 42077270.

Data Availability Statement

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 conflict of interest.

References

  1. Shi, G.H. Discontinuous Deformation Analysis—A New Numerical Model for the Statics and Dynamics of Block Systems. Ph.D. Thesis, Department of Civil Engineering University of California, Berkeley, CA, USA, 1988. [Google Scholar]
  2. MacLaughlin, M.M.; Doolin, D.M. Review of validation of the discontinuous deformation analysis (DDA) method. Int. J. Numer. Anal. Meth. Geomech. 2006, 30, 271–305. [Google Scholar] [CrossRef]
  3. Jing, L.R. Formulation of discontinuous deformation analysis (DDA)—an implicit discrete element model for block systems. Eng. Geol. 1998, 49, 371–381. [Google Scholar] [CrossRef]
  4. Yang, Y.T.; Xu, D.D.; Liu, F.; Zheng, H. Modeling the entire progressive failure process of rock slopes using a strength-based criterion. Comput. Geotech. 2020, 126, 103726. [Google Scholar] [CrossRef]
  5. Chen, K.T.; Wu, J.H. Simulating the failure process of the Xin Mo landslide using discontinuous deformation analysis. Eng. Geol. 2018, 239, 269–281. [Google Scholar] [CrossRef]
  6. Hatzor, Y.H.; Arzi, A.A.; Zaslavsky, Y.; Shapira, A. Dynamic stability analysis of jointed rock slopes using the DDA method: King Herod’s Palace, Masada, Israel. Int. J. Rock Mech. Min. Sci. 2004, 41, 813–832. [Google Scholar] [CrossRef]
  7. Fu, X.D.; Sheng, Q.; Zhang, Y.H.; Chen, J.; Zhang, S.K.; Zhang, Z.P. Computation of the safety factor for slope stability using discontinuous deformation analysis and the vector sum method. Comput. Geotech. 2017, 92, 68–76. [Google Scholar] [CrossRef]
  8. Ning, Y.J.; Yang, J.; An, X.M.; Ma, G.W. Modelling rock fracturing and blast-induced rock mass failure via advanced discretization within the discontinuous deformation analysis framework. Comput. Geotech. 2011, 38, 40–49. [Google Scholar] [CrossRef]
  9. Tang, C.A.; Tang, S.B.; Gong, B.; Bai, H.M. Discontinuous deformation and displacement analysis: From continuous to discontinuous. Sci. China Tech. Sci. 2015, 58, 1567–1574. [Google Scholar] [CrossRef]
  10. Chen, Y.J.; Zhang, X.; Zhu, W.S.; Wang, W. Modified discontinuous deformation analysis for rock failure: Crack propagation. Geomech. Eng. 2018, 14, 325–336. [Google Scholar]
  11. Liu, X.Y.; Zhao, Z.Y.; Ma, S.Q. Simulating stress wave with flat-top partition of unity based high-order discontinuous deformation analysis. Eng. Anal. Bound. Elem. 2018, 91, 110–123. [Google Scholar] [CrossRef]
  12. Chen, G.Q.; He, M.C.; Fan, F.S. Rock burst analysis using DDA numerical simulation. Int. J. Geomech. 2018, 18, 04018001. [Google Scholar] [CrossRef]
  13. Ning, Y.; Yang, Z.; Wei, B.; Gu, B. Advances in two-dimensional discontinuous deformation analysis for rock-mass dynamics. Int. J. Geomech. 2016, 17, E6016001. [Google Scholar] [CrossRef]
  14. Cai, Y.; Liang, G.; Shi, G.H. Studying impact problem by LDDA method. In Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media; TSI Press: San Antonio, TX, USA, 1996. [Google Scholar]
  15. Bao, H.R.; Zhao, Z.Y.; Tian, Q. On the implementation of augmented Lagrangian method in the two-dimensional discontinuous deformation analysis. Int. J. Numer. Anal. Meth. Geomech. 2014, 38, 551–571. [Google Scholar] [CrossRef]
  16. Zheng, H.; Jiang, W. Discontinuous deformation analysis based on complementary theory. Sci. China Ser. E Technol. Sci. 2009, 52, 2547–2554. [Google Scholar] [CrossRef]
  17. Jiang, W.; Zheng, H. Discontinuous deformation analysis based on variational inequality theory. Int. J. Comput. Methods 2011, 8, 193–208. [Google Scholar] [CrossRef]
  18. Lin, S.; Xie, Z. Performance of DDA time integration. Sci. China Technol. Sci. 2015, 58, 1558–1566. [Google Scholar] [CrossRef]
  19. Jiao, Y.Y.; Zhang, X.L.; Zhao, J.; Liu, Q.S. Viscous boundary of DDA for modeling stress wave propagation in jointed rock. Int. J. Rock Mech. Min. Sci. 2007, 44, 1070–1076. [Google Scholar] [CrossRef]
  20. Jiang, Q.; Chen, Y.; Zhou, C.; Yeung, M.-C.R. Kinetic Energy Dissipation and Convergence Criterion of Discontinuous Deformations Analysis (DDA) for Geotechnical Engineering. Rock Mech. Rock Eng. 2012, 46, 1443–1460. [Google Scholar] [CrossRef]
  21. Beyabanaki, S.A.R.; Jafari, A.; Yeung, M.R. High-order three-dimensional discontinuous deformation analysis (3-D DDA). Int. J. Numer. Meth. Biomed. Eng. 2010, 26, 1522–1547. [Google Scholar] [CrossRef]
  22. Grayeli, R.; Hatami, K. Implementation of the finite element method in the three-dimensional discontinuous deformation analysis (3D-DDA). Int. J. Numer. Anal. Meth. Geomech. 2008, 32, 1883–1902. [Google Scholar] [CrossRef]
  23. Wang, J.Q.; Lin, G. Static and dynamic stability analysis using 3D-DDA with incision body scheme. Earthq. Eng. Eng. Vib. 2006, 5, 273–284. [Google Scholar] [CrossRef]
  24. Liu, J.; Geng, Q.D.; Kong, X.J. A fast common plane identification algorithm for 3D contact problems. Comput. Geotech. 2009, 36, 41–51. [Google Scholar] [CrossRef]
  25. Ali, R.K.; Ahmad, J.; Wu, J.H. A new algorithm to identify contact patterns between convex blocks for three-dimensional discontinuous deformation analysis. Comput. Geotech. 2008, 35, 746–759. [Google Scholar]
  26. Zhu, H.; Wu, W.; Chen, J.; Ma, G.; Liu, X.; Zhuang, X. Integration of three dimensional discontinuous deformation analysis (DDA) with binocular photogrammetry for stability analysis of tunnels in blocky rock mass. Tunn. Undergr. Space Technol. 2016, 51, 30–40. [Google Scholar] [CrossRef]
  27. Jiao, Y.Y.; Huang, G.H.; Zhao, Z.Y.; Zheng, F.; Wang, L. An improved three-dimensional spherical DDA model for simulating rock failure. Sci. China Tech. Sci. 2015, 58, 1533–1541. [Google Scholar] [CrossRef]
  28. Shi, G.H. Contact theory. Sci. China Tech. Sci. 2015, 58, 1450–1496. [Google Scholar] [CrossRef]
  29. Song, X.Y.; Huang, D.; Zeng, B. GPU-based parallel computation for discontinuous deformation analysis (DDA) method and its application to modelling earthquake-induced landslide. Comput. Geotech. 2017, 86, 80–94. [Google Scholar] [CrossRef]
  30. Yang, Y.; Xu, D.; Zheng, H. Explicit Discontinuous Deformation Analysis Method with Lumped Mass Matrix for Highly Discrete Block System. Int. J. Geomech. 2018, 18, 04018098. [Google Scholar] [CrossRef]
  31. MacLaughlin, M.M.; Sitar, N. Rigid body rotations in DDA. In Proceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, Berkeley, CA, USA, 12–14 June 1996; TSI Press: San Antonio, TX, USA, 1996; pp. 620–635. [Google Scholar]
  32. Cheng, Y.M.; Zhang, Y.H. Rigid body rotation and block internal discretization in DDA analysis. Int. J. Numer. Anal. Meth. Geomech. 2000, 24, 567–578. [Google Scholar] [CrossRef]
  33. Jiang, W.; Zheng, H. An efficient remedy for the false volume expansion of DDA when simulating large rotation. Comput. Geotech. 2015, 70, 18–23. [Google Scholar] [CrossRef]
  34. Koo, C.Y.; Chern, J.C. Modification of the DDA method for rigid block problems. Int. J. Rock Mech. Mining. Sci. 1998, 35, 683–693. [Google Scholar] [CrossRef]
  35. Ke, T.C. The Issue of Rigid-Body Rotation in DDA. In Proceedings of the First International Forum on Discontinuous Deformation Analysis (DDA) and Simulations of Discontinuous Media, Berkeley, CA, USA, 12–14 June 1996; TSI Press: San Antonio, TX, USA, 1996; pp. 318–326. [Google Scholar]
  36. Jiang, W.; Wang, Y.H.; Sun, G.H.; Song, P.C.; Mao, C. Discontinuous Deformation Analysis Based on Wachspress Interpolation Function for Detailed Stress Distribution. Math. Probl. Eng. 2018, 2018, 4079035. [Google Scholar] [CrossRef]
  37. Dasgupta, G. Interpolants within convex polygon: Wachspress′ shape functions. J. Aerosp. Eng. 2003, 16, 1–8. [Google Scholar] [CrossRef] [Green Version]
  38. Sukumar, N. Construction of polygonal interpolants: A maximum entropy interpolation. Int. J. Numer. Meth. Eng. 2004, 61, 2159–2181. [Google Scholar] [CrossRef]
  39. Jiang, W.; Xu, J.J.; Zheng, H.; Wang, Y.H.; Sun, G.H.; Yang, Y.T. Novel displacement function for discontinuous deformation analysis based on mean value coordinates. Int. J. Numer. Meth. Eng. 2020, 121, 4768–4792. [Google Scholar] [CrossRef]
  40. Joshi, P.; Meyer, M.; DeRose, T.; Green, B.; Sanocki, T. Harmonic coordinates for character articulation. ACM Trans. Graph. 2007, 26, 1–9. [Google Scholar] [CrossRef]
  41. Talischi, C.; Pereira, A.; Menezes, I.F.M.; Paulino, G.H. Gradient correction for polygonal and polyhedral finite elements. Int. J. Numer. Methods Eng. 2015, 102, 728–747. [Google Scholar] [CrossRef]
  42. Talischi, C.; Paulino, G.H. Addressing integration error for polygonal finite elements through polynomial projections: A patch test connection. Math. Models Methods Appl. Sci. 2014, 24, 1701–1727. [Google Scholar] [CrossRef]
  43. Beirão da Veiga, L.; Brezzi, F.; Cangiani, A.; Manzini, G.; Marini, L.D.; Russo, A. Basic principles of virtual element methods. Math. Models Methods Appl. Sci. 2013, 23, 199–214. [Google Scholar] [CrossRef]
  44. Ahmad, B.; Alsaedi, A.; Brezzi, F.; Marini, L.D.; Russo, A. Equivalent projectors for virtual element methods. Comput. Math. Appl. 2013, 66, 376–391. [Google Scholar] [CrossRef]
  45. da Veiga LBo Brezzi, F.; Marini, L.D. Virtual Elements for Linear Elasticity Problems. SIAM J. Numer. Anal. 2013, 51, 794–812. [Google Scholar] [CrossRef] [Green Version]
  46. Beirão da Veiga, L.; Brezzi, F.; Marini, L.D.; Russo, A. The hitchhiker′s guide to the virtual element method. Math. Models Methods Appl. Sci. 2014, 24, 1541–1573. [Google Scholar] [CrossRef] [Green Version]
  47. Park, K.; Chi, H.; Paulino, G.H. On nonconvex meshes for elastodynamics using virtual element methods with explicit time integration. Comput. Methods Appl. Mech. Eng. 2019, 356, 669–684. [Google Scholar] [CrossRef]
  48. Gain, A.L.; Talischi, C.; Paulino, G.H. On the Virtual Element Method for three-dimensional linear elasticity problems on arbitrary polyhedral meshes. Comput. Methods Appl. Mech. Eng. 2014, 282, 132–160. [Google Scholar] [CrossRef] [Green Version]
  49. Wriggers, P.; Rust, W.T.; Reddy, B.D. A virtual element method for contact. Comput. Mech. 2016, 58, 1039–1050. [Google Scholar] [CrossRef]
  50. Beirão da Veiga, L.; Brezzi, F.; Dassi, F.; Marini, L.D.; Russo, A. Virtual Element approximation of 2D magnetostatic problems. Comput. Methods Appl. Mech. Eng. 2017, 327, 173–195. [Google Scholar] [CrossRef]
  51. Wriggers, P.; Reddy, B.D.; Rust, W.; Hudobivnik, B. Efficient virtual element formulations for compressible and incompressible finite deformations. Comput. Mech. 2017, 60, 253–268. [Google Scholar] [CrossRef]
  52. Chi, H.; da Veiga, L.B.; Paulino, G. Some basic formulations of the virtual element method (VEM) for finite deformations. Comput. Methods Appl. Mech. Eng. 2017, 318, 148–192. [Google Scholar] [CrossRef]
Figure 1. A typical block system and an individual block in Discontinuous Deformation Analysis (DDA).
Figure 1. A typical block system and an individual block in Discontinuous Deformation Analysis (DDA).
Materials 14 01252 g001
Figure 2. A polygon and its boundary division.
Figure 2. A polygon and its boundary division.
Materials 14 01252 g002
Figure 3. Triangles involved in the calculation of Ni using Wachspress interpolation.
Figure 3. Triangles involved in the calculation of Ni using Wachspress interpolation.
Materials 14 01252 g003
Figure 4. The notation li and ni in H.
Figure 4. The notation li and ni in H.
Materials 14 01252 g004
Figure 5. A typical contact pair at (a) the beginning and (b) the end of the current time step.
Figure 5. A typical contact pair at (a) the beginning and (b) the end of the current time step.
Materials 14 01252 g005
Figure 6. (a) A triangular block on the foundation and (b) its movement under a horizontal point load.
Figure 6. (a) A triangular block on the foundation and (b) its movement under a horizontal point load.
Materials 14 01252 g006
Figure 7. (a) The configuration of the sliding problem and (b) the relative displacement errors of the proposed method and DDA with the post-adjustment strategy.
Figure 7. (a) The configuration of the sliding problem and (b) the relative displacement errors of the proposed method and DDA with the post-adjustment strategy.
Materials 14 01252 g007
Figure 8. (a) The variations in minimum principal stress direction with the friction angles and (b) the difference in minimum principal stress direction of the two methods.
Figure 8. (a) The variations in minimum principal stress direction with the friction angles and (b) the difference in minimum principal stress direction of the two methods.
Materials 14 01252 g008
Figure 9. A surrounding rock model.
Figure 9. A surrounding rock model.
Materials 14 01252 g009
Figure 10. Configurations and main stresses after 100 time steps under two different friction angles: (a) φ = 30° and (b) φ = 0°.
Figure 10. Configurations and main stresses after 100 time steps under two different friction angles: (a) φ = 30° and (b) φ = 0°.
Materials 14 01252 g010
Figure 11. The block wall model.
Figure 11. The block wall model.
Materials 14 01252 g011
Figure 12. The configurations and main stresses after 7 s using the proposed method.
Figure 12. The configurations and main stresses after 7 s using the proposed method.
Materials 14 01252 g012
Figure 13. The configurations and main stress vectors after 7 s by DDA with the post-adjustment strategy.
Figure 13. The configurations and main stress vectors after 7 s by DDA with the post-adjustment strategy.
Materials 14 01252 g013
Table 1. The penetrations of the contact pair point 7–edge 14 by DDA with the post-adjustment strategy and the new method.
Table 1. The penetrations of the contact pair point 7–edge 14 by DDA with the post-adjustment strategy and the new method.
Time StepThe Original DDA with the Post-Adjustment StrategyThe New MethodThe Rotational Angle r0
Penetrations in the Equilibrium Equations (m)Penetrations in the Updated Configurations (m)Penetrations in the Equilibrium Equations (m)Penetrations in the Updated Configurations (m)
18.3289 × 1098.3288 × 1098.32890 × 1098.32890 × 109−0.000001
22.5026 × 1082.5026 × 1082.50260 × 1082.50260 × 108−0.000000
37.4918 × 1087.4918 × 1087.49180 × 1087.49180 × 108−0.000000
43.3772 × 1071.4587 × 1073.37720 × 1073.37720 × 107−0.001073
51.0094 × 106−7.2045 × 1071.00970 × 1061.00970 × 106−0.003218
63.0067 × 106−1.8379 × 1063.01240 × 1063.01240 × 106−0.005375
78.9005 × 106−7.2391 × 1078.92310 × 1068.92310 × 106−0.007554
82.6114 × 1059.8896 × 1062.62290 × 1052.62290 × 105−0.009768
97.5032 × 1055.0108 × 1057.54400 × 1057.54400 × 105−0.012047
102.0631 × 1041.7005 × 1042.07880 × 1042.07880 × 104−0.014444
114.9419 × 1044.4322 × 1044.97090 × 1044.97090 × 104−0.017006
127.6165 × 1046.9448 × 1047.63320 × 1047.63320 × 104−0.019372
136.5260 × 1045.7022 × 1046.50170 × 1046.50170 × 104−0.021275
143.2463 × 1042.2251 × 1043.22710 × 1043.22710 × 104−0.023478
152.6022 × 1041.2772 × 1042.68970 × 1042.68970 × 104−0.026494
165.7111 × 1043.9996 × 1045.96770 × 1045.96770 × 104−0.029810
179.6448 × 1047.5227 × 1049.91470 × 1049.91470 × 104−0.032847
189.2970 × 1046.7950 × 1049.13440 × 1049.13440 × 104−0.035286
193.3880 × 1044.7618 × 1052.94330 × 1042.94330 × 104−0.037664
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Luo, H.; Sun, G.; Liu, L.; Jiang, W. Vertex Displacement-Based Discontinuous Deformation Analysis Using Virtual Element Method. Materials 2021, 14, 1252. https://0-doi-org.brum.beds.ac.uk/10.3390/ma14051252

AMA Style

Luo H, Sun G, Liu L, Jiang W. Vertex Displacement-Based Discontinuous Deformation Analysis Using Virtual Element Method. Materials. 2021; 14(5):1252. https://0-doi-org.brum.beds.ac.uk/10.3390/ma14051252

Chicago/Turabian Style

Luo, Hongming, Guanhua Sun, Lipeng Liu, and Wei Jiang. 2021. "Vertex Displacement-Based Discontinuous Deformation Analysis Using Virtual Element Method" Materials 14, no. 5: 1252. https://0-doi-org.brum.beds.ac.uk/10.3390/ma14051252

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop