# Discontinuous Deformation Analysis Coupling with Discontinuous Galerkin Finite Element Methods for Contact Simulations

Authors: Yue Sun; Xiangchu Feng; Jun Xiao; Ying Wang
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

A novel coupling scheme is presented to combine the discontinuous deformation analysis (DDA) and the interior penalty Galerkin (IPG) method for the modeling of contacts. The simultaneous equilibrium equations are assembled in a mixed strategy, where the entries are derived from both discontinuous Galerkin variational formulations and the strain energies of DDA contact springs. The contact algorithms of the DDA are generalized for element contacts, including contact detection criteria, open-close iteration, and contact submatrices. Three representative numerical examples on contact problems are conducted. Comparative investigations on the results obtained by our coupling scheme, ANSYS, and analytical theories demonstrate the accuracy and effectiveness of the proposed method.


---

## Body

## 1. Introduction

In recent years, the discontinuous deformation analysis (DDA) [1, 2] has become increasingly researched in a vast variety of geotechnical engineering applications ranging from seismic landslides [3, 4] to hydrofracturing for shale gas exploitation [5–8]. At either end of this spectrum, reliable and robust modeling of contact interactions have vital significances in deciding the initiation of crack tips as well as their growths. So far, extensive endeavors have been dedicated to boost the efficiency of DDA contact algorithms [9–18]. However, since block stresses are assumed as constants, the contact pressures figured out by this method are usually of low accuracy. For this reason, a multitude of approaches have been proposed to facilitate DDA block deformability, such as introducing artificial joints [19, 20], adopting higher-order polynomial basis functions [21], incorporating finite element meshes [22–26], and coupling the DDA and finite element methods (FEMs) [27–30].Compared with the other enriched DDA approaches, coupling methods can take the advantages of both the DDA for contacts and the FEM [31, 32] for accurate block deformations. There have existed three kinds of coupling schemes at present: (a) Zhang et al. [28] applied the DDA and the FEM on two different blocks, respectively, and introduced a line block to measure penetration depths; (b) the authors of literatures [27, 29] alternatively applied the DDA and the FEM to figure out contact forces and block deformations; (c) Zhao and Gu [30] used a stress recovery procedure to refine DDA results. The limitation of the scheme (a) lies in its very strict requirement on block displacements in one step, and this method can only be used in static problems. Since the scheme (b) needs to solve both DDA and FEM equations and transfer their results to each other repeatedly in every time step, the additional FEM calculation leads to higher computational cost and hence is detrimental to its efficiency. The scheme (c) also has the disadvantage similar as the scheme (b) due to the stress recovery phase, in which nodal force balance equations have to be solved once a DDA computation finishes.The objective of this article is to present a novel coupling scheme that attempts to combine the DDA and discontinuous Galerkin (DG) FEMs [33, 34]. DG FEMs have plenty of outstanding features, such as supporting irregularly shaped elements and nonconforming meshes, allowing discontinuities in both trial and test function spaces, local conservations of momentums and masses, and well established mathematical theories [35]. Particularly, the DG FEM employed in this paper is known as the interior penalty Galerkin method (IPG) [36–38]. For clarity, we refer to our scheme hereafter as the IPG-DDA.In the proposed scheme, the DG-FEM is implemented on the blocks whose stresses are to be refined by subdividing them into triangular elements, while each of the rest intact blocks is regarded as a block-type element. Consequently, the original multiblock system becomes assembles of elements. The contacts between blocks transfer to the ones along the element edges which compose the boundaries of the original blocks. Contact springs are applied at every contact point and controlled by the open-close iteration similar as in the standard DDA. During each time step, the global equilibrium equations are assembled in a mixed strategy: the DG-FEM contributes the submatrices derived from the Galerkin-based variational formulations in broken Sobolev spaces, while the DDA gives the submatrices generated by minimizing the strain energies of the contact springs between exterior element edges. In this way, refined block deformations can be obtained by solving the global equations, and contact forces need not be figured out alone. Two auxiliary criteria are also proposed to determine the types of element contacts for the IPG-DDA.The differences between our scheme and the aforementioned schemes (a), (b), and (c) lie in the coupling methodology and the element types adopted in them. For the first difference, the IPG-DDA combines the DDA and the FEM in an integrated framework but does not apply them alternatively or, respectively, on different block domains. Therefore, we avoid solving additional FEM equations and transferring the results forth and back to the DDA computation. For the second difference, the IPG method is a kind of stable mixed FEMs rather than the conforming FEM adopted in the existed coupling schemes. The advantages of the IPG method include accurate, no-locking, low computing cost, and the great flexibility in meshing. Plenty of numerical experiments on contact issues are conducted to verify the presented IPG-DDA method. The results are compared with analytical theories and ANSYS solutions, and our coupling scheme is demonstrated to be accurate and effective.Although the methods known as the nodal-based DDA [22–26] are also regarded as an another sort of coupling schemes with the FEM in some literatures, there exist important differences between our scheme and these methods. Nodal-based DDA incorporates finite element type meshes on blocks and uses nodal displacements as the unknowns to be determined. Hence, their global stiffness matrices require the elements having the same degrees of freedom (DOFs) to avoid rank-deficiency. In contrast to the nodal-based DDA, our method permits the elements with any shapes, even original DDA blocks. Additionally, nodal stresses are assumed to be continuous in the nodal-based DDA. However, so strong requirement restricts element behaviors from large strains and distortions. Our scheme allows completely discontinuous displacements and stresses across element edges, thus free from those restrictions.


## 2. Linear Elasticity Dynamics on a Single Block

Consider the block system containingN linearly elastic polygonal blocks. In this section, we give the variational formulations of the IPG method on single block domains. In the next section, the submatrices derived from these variational formulations are assembled to the global equilibrium equations.


### 2.1. Problem Definition

Assume the blockΩ with the boundary Γ has the mass density ρ. Ω is subjected to the body force b per unit area. Γ consists of two mutually disjoint parts Γu and Γt, where the prescribed displacement u¯ and the surface traction t¯ are defined, respectively. Use u to denote the displacement and σ for the Cauchy stress tensor of the block Ω; then the dynamics of the block Ω at the time t∈0,∞ can be stated as(1)ρu¨=∇·σ+b,inΩ,(2)u=u¯,onΓu,(3)σn=t¯,onΓt,where ∇ is the gradient operator and n represents the unit outward normal vector of the boundary Γ. Use ε to indicate the strain tensor, and there exist the following relationships:(4)σ=C:ε,inΩ,(5)ε=12∇u+12∇uT,inΩ,where C is a fourth-order elastic tensor, which is symmetric, positive definite, and point-wisely constant in Ω.


### 2.2. Finite Element Discretization

The IPG method discretizes a block into elements, then approximates primary variables by discontinuous shape functions, and weakly enforces interelement continuity and essential boundary constraints by penalty techniques. LetTh={E1,E2,…,EΛ} be a nondegenerate (the nondegeneracy means that there exists ri>0 such that if hi=diamEi then Ei contains a ball of radius rihi in its interior. The operator diam(Ei) is the supremum of the distance between any two points in Ei) triangular subdivision of the block Ω. Denote the edges of all elements in Th by E=e1,e2,…,ePα,ePα+1,…,eQα, where e1,e2,…,ePα⊂Ω and Γ=ePα+1,…,eQα. As shown in Figure 1, we call eξ as an interior edge if 1⩽ξ⩽Pα or an exterior edge if Pα+1⩽ξ⩽Qα. Denote the edges of the element E by ∂E. Then the finite element subspace on Th can be defined as [34](6)DkTh=v∈L2Ω:vE∈PkE,∀E∈Th,where L2Ω denotes the space of square-integrable functions and PkE is the space of the polynomials of the degree less than or equal to the positive integer k. DkTh is full of discontinuous piecewise polynomials and is used as the test function space below.Figure 1
Surface-surface contact between two blocks.


### 2.3. DG-FEM Weak Form

According to the IPG method [34, 36, 37], the variational formulation of problem (1)–(3) is to find u∈DkTh, such that(7)au,v=Lv,∀v∈DkTh,where the bilinear form au,v:DkTh×DkTh→R:(8)au,v=∑E∈Th∫Eρu¨·vdE+∑E∈Th∫Eσu:εvdE-∑e∈E∖Γ∪Γu∫eσune·vdl-∑e∈E∖Γ∪Γu∫eσvne·udl+∑e∈E∖Γ∪ΓuδpGe∫eu·vdl,and the linear form Lv:DkTh→R:(9)Lv=∑E∈Th∫Eb·vdE+∑e∈Γt∫et¯·vdl+∑e∈Γu∫eσvne·u¯dl+∑e∈ΓuδpGe∫eu¯·vdl,where δp and G denote the penalty parameter and the shear modulus, respectively, e is the length of the element edge e∈E, and ne denotes the unit outward normal vector to e; the average operator · and the jump operator · across e are defined as(10)w=12wE+e+wE-e,e∈E∖Γwie,e∈Γ,(11)w=wE+e-wE-e,e∈E∖Γwie,e∈Γ,where e is the common edge of the elements E+ and E-.Substitute (8) and (9) into (7), and we have(12)∑E∈Thaρ+∑E∈Thaσ+∑e∈Sas=∑E∈Thlb+∑e∈Γtlt+∑e∈Γulu,where S≡E∖Γ∪Γu, and(13)aρ=∫E∈Thρu¨·vdE,(14)aσ=∫E∈Thσu:εvdE,(15)as=-∫e∈Sσune·vdl-∫e∈Sσvne·udl+δpGe∫e∈Su·vdl,(16)lb=∫E∈Thb·vdE,(17)lt=∫e∈Γtt¯·vdl,(18)lu=∫e∈Γuσvne·u¯dl+δpGe∫e∈Γuu¯·vdl,where the operator “·” indicates the dot product of two vectors.


## 2.1. Problem Definition

Assume the blockΩ with the boundary Γ has the mass density ρ. Ω is subjected to the body force b per unit area. Γ consists of two mutually disjoint parts Γu and Γt, where the prescribed displacement u¯ and the surface traction t¯ are defined, respectively. Use u to denote the displacement and σ for the Cauchy stress tensor of the block Ω; then the dynamics of the block Ω at the time t∈0,∞ can be stated as(1)ρu¨=∇·σ+b,inΩ,(2)u=u¯,onΓu,(3)σn=t¯,onΓt,where ∇ is the gradient operator and n represents the unit outward normal vector of the boundary Γ. Use ε to indicate the strain tensor, and there exist the following relationships:(4)σ=C:ε,inΩ,(5)ε=12∇u+12∇uT,inΩ,where C is a fourth-order elastic tensor, which is symmetric, positive definite, and point-wisely constant in Ω.


## 2.2. Finite Element Discretization

The IPG method discretizes a block into elements, then approximates primary variables by discontinuous shape functions, and weakly enforces interelement continuity and essential boundary constraints by penalty techniques. LetTh={E1,E2,…,EΛ} be a nondegenerate (the nondegeneracy means that there exists ri>0 such that if hi=diamEi then Ei contains a ball of radius rihi in its interior. The operator diam(Ei) is the supremum of the distance between any two points in Ei) triangular subdivision of the block Ω. Denote the edges of all elements in Th by E=e1,e2,…,ePα,ePα+1,…,eQα, where e1,e2,…,ePα⊂Ω and Γ=ePα+1,…,eQα. As shown in Figure 1, we call eξ as an interior edge if 1⩽ξ⩽Pα or an exterior edge if Pα+1⩽ξ⩽Qα. Denote the edges of the element E by ∂E. Then the finite element subspace on Th can be defined as [34](6)DkTh=v∈L2Ω:vE∈PkE,∀E∈Th,where L2Ω denotes the space of square-integrable functions and PkE is the space of the polynomials of the degree less than or equal to the positive integer k. DkTh is full of discontinuous piecewise polynomials and is used as the test function space below.Figure 1
Surface-surface contact between two blocks.


## 2.3. DG-FEM Weak Form

According to the IPG method [34, 36, 37], the variational formulation of problem (1)–(3) is to find u∈DkTh, such that(7)au,v=Lv,∀v∈DkTh,where the bilinear form au,v:DkTh×DkTh→R:(8)au,v=∑E∈Th∫Eρu¨·vdE+∑E∈Th∫Eσu:εvdE-∑e∈E∖Γ∪Γu∫eσune·vdl-∑e∈E∖Γ∪Γu∫eσvne·udl+∑e∈E∖Γ∪ΓuδpGe∫eu·vdl,and the linear form Lv:DkTh→R:(9)Lv=∑E∈Th∫Eb·vdE+∑e∈Γt∫et¯·vdl+∑e∈Γu∫eσvne·u¯dl+∑e∈ΓuδpGe∫eu¯·vdl,where δp and G denote the penalty parameter and the shear modulus, respectively, e is the length of the element edge e∈E, and ne denotes the unit outward normal vector to e; the average operator · and the jump operator · across e are defined as(10)w=12wE+e+wE-e,e∈E∖Γwie,e∈Γ,(11)w=wE+e-wE-e,e∈E∖Γwie,e∈Γ,where e is the common edge of the elements E+ and E-.Substitute (8) and (9) into (7), and we have(12)∑E∈Thaρ+∑E∈Thaσ+∑e∈Sas=∑E∈Thlb+∑e∈Γtlt+∑e∈Γulu,where S≡E∖Γ∪Γu, and(13)aρ=∫E∈Thρu¨·vdE,(14)aσ=∫E∈Thσu:εvdE,(15)as=-∫e∈Sσune·vdl-∫e∈Sσvne·udl+δpGe∫e∈Su·vdl,(16)lb=∫E∈Thb·vdE,(17)lt=∫e∈Γtt¯·vdl,(18)lu=∫e∈Γuσvne·u¯dl+δpGe∫e∈Γuu¯·vdl,where the operator “·” indicates the dot product of two vectors.


## 3. Global Computing Scheme

Consider theN-block system Ω1,Ω2,…,ΩN, Ωα⊂R2. Subdivide the blocks whose stresses are to be refined, and the rest intact blocks are regarded as block-type elements. Assume M elements are created from the considering block system, including E1, E2,…,EM. Numerically discretize the displacement ui of the element Ei by(19)ui=uix,yvix,y=Tix,yDi,∀x,y∈Ei,where Di denotes the unknown vector variable to be determined on the element Ei. We adopt the following Tix,y and Di:(20)Tix,y=10y~i-yx-x~i0y-y~i201x-x~i0y-y~ix-x~i2,(21)Di=u0v0r0εxεyγxyT,where x~i,y~i denotes the barycenter coordinate of the element Ei, u0, v0 are the translational displacements at x~,y~, r0 is the rotation at x~i,y~i, and εx, εy, γxy denote the strain components of the element Ei.In every time step, the following global equilibrium equations are formed in a mixed strategy:(22)K11K12⋯K1MK21K22⋯K2M⋮⋮⋱⋮KM1KM2⋯KMMD1D2⋮DM=F1F2⋮FM,where the stiffness matrix is partitioned by the submatrices Kij∈R6×6, the right hand side consists of Fi∈R6×1, and Di∈R6×1 contains the unknowns to be determined on the ith element, i,j=1,2,…,M. The mixed strategy for assembling (22) can be explained through the example as shown in Figure 2. In Figure 2(a), there exist three contact pairs, including (Block 1, Block 2), (Block 1, Block 3), and (Block 2, Block 3). In this case, the stiffness matrix of (22) has the block structure as illustrated in Figure 2(b), where shadowed diagonal partitions are derived from both the variational formulation equation (12) and the minimizations of DDA contact spring strain energies, while the off-diagonal partitions are contributed by only contact springs.Figure 2
The structure of the stiffness matrix of IPG-DDA global simultaneous equilibrium equations in the case of three blocks.
(a)
Stacked three blocks and the implemented meshes (b)
Structure of the stiffness matrixIt is demonstrated that the coupling begins in the phase of assembling the global equilibrium equations in our method. Such a mixed strategy is very different from the existing coupling schemes in [27–30]. The latter ones need to solve DDA and FEM equations alternatively, and in time-dependent problems they also have to transfer the solutions of the two methods with each other repeatedly. The IPG-DDA simultaneously figures out block deformations and contact displacements through the unified global equilibrium equations (22).The entriesKij and Fi in (22) are derived below. Sections 3.1 to 3.6 are based on the DG-FEM variational formulation (12). Sections 3.7 and 3.8 generalize DDA contact submatrices [1].


### 3.1. Inertia Force

In the time step intervalt,t+Δt, substituting (19) into (13) leads to(23)aρ=DiT∫EiρTiTTidE∂2∂t2Di≡DiTMi∂2∂t2Di,where the area integral in the parentheses can be calculated by the simplex integral formulations [39]. With the help of Newmark time integral scheme [40], we can solve the second derivative of Di in (23) as(24)aρ=1βΔt2DiTMiDi-1βΔtDiTMiVi0-12β-1DiTMiAi0,where γ and β indicate Newmark parameters and Vi0 and Ai0 denote the initial velocity and acceleration of the current time step. Notice that the initial velocity and acceleration in one time step are actually the finial velocity and acceleration in the last time step. Therefore we can use the following recursive formulations to update Vi0 and Ai0 in every time step:(25)Ai1=1βΔt2Di-1βΔtVi0-12β-1Ai0,Vi1=γβΔtDi+1-γβVi0+1-γ2βΔtAi0,where Vi1 and Ai1 refer to the finial velocity and acceleration, respectively. In the numerical examples of this article, we choose γ=1 and β=0.5 to avoid numerical damping and to ensure this time integration method to be unconditionally stable.Since the IPG method is based on the Galerkin weighted residual method, the global equilibrium equations have to hold for all test functionsvi, and therefore the vector DiT in (24) can be divided out from the equations. So we assemble 1/βΔt2Mi on the right side of (24) to [Kii] in the stiffness matrix of the global equilibrium equations (22), while -1/βΔtMiVi0 and -1/2β-1MiAi0 are added to Fi of the right hand side of (22). The submatrices derived below are assembled to the global equilibrium equations due to the same reason.


### 3.2. Element Stiffness

Denote the vector forms of the stress tensor and the strain tensor on the elementEi by σi and εi, respectively. Substituting (21) into (4) and (5), respectively, leads to(26)εi=εxεyγxy=000100000010000001Di≡BDi,σi=σxσyτxy=E¯1-ν¯21ν¯0ν¯10001-ν¯2εi≡CiBDi,where E¯=E and ν¯=ν for plane stress problems and E¯=E/(1-ν2) and ν¯=ν/(1-ν) for plane strain problems, in which the constants E and ν indicate Young’s modulus and Poisson’s ratio of the element Ei, respectively. The matrix [Ci] above is known as the material matrix of Ei. For the three-dimensional case rather than the two-dimensional one as concerned in this paper, contact points, edges, and planes have to be taken into account and detailed in different manners.Then substitute (19) and (26) into (14), and we have(27)aσ=∫EiDiTBTCiBDidEi=DiTEiBTCiBDi,where Ei is the area of the element Ei. The matrix EiBTCiB is assembled to the 6×6 submatrix Kii in the stiffness matrix of the global equilibrium equations (22).


### 3.3. Body Force

Assume the elementEi bearing the body force bi=bx,byT. Substituting (19) into (16) gives(28)lb=DiT∫EiTTbxbydEi,where the integral on the right of the above equality is assembled to the right hand side Fi of the global simultaneous equilibrium equations (22).


### 3.4. Interior Interface Stiffness

Assume the interior edgee as the common boundary of the elements Ei and Ej. Rewrite as in (15) into the summation below:(29)as=as,1+as,2+as,3,where(30)as,1=-∫eσune·vdl,as,2=-∫eσvne·udl,as,3=δpGe∫eu·vdl.Using the notations of average and jump operators defined in (10) and (11), respectively, we can expand the expressions above into the formulations below:(31)as,1=-12∫eσuine·vidl+12∫eσujne·vidl+12∫eσuine·vjdl-12∫eσujne·vjdl,(32)as,2=-12∫eσvine·uidl+12∫eσvjne·uidl+12∫eσvine·ujdl-12∫eσvjne·ujdl,(33)as,3=δpGe∫eui·vidl-δpGe∫eui·vjdl-δpGe∫euj·vidl+δpGe∫euj·vjdl.Then substitute (19) into (31), (32), and (33), respectively, and we have(34)as,1=-12DiT∫eTiTdlNeCBDi+12DjT∫eTiTdlNeCBDj+12DiT∫eTjTdlNeCBDi-12DjT∫eTjTdlNeCBDj,as,2=-12DiTBTCTNeT∫eTiTdlDi+12DjTBTCTNeT∫eTidlDi+12DiTBTCTNeT∫eTjTdlDj-12DjTBTCTNeT∫eTjTdlDj,as,3=δpGeDiT∫eTiTTidlDi-δpGeDjT∫eTjTTidlDi-δpGeDiT∫eTiTTjdlDj+δpGeDjT∫eTjTTjdlDj,where the matrix Ne is spanned by the unit outward normal vector nx,ny to the edge e:(35)Ne=nx0ny0nynx.Actually, Ne transfers the surface traction f¯x,f¯yT to the combination of stress components in the following way:(36)f¯xf¯y=σne=σxτxyτxyσynxny=Neσxσyτxy,where σ denotes the stress tensor referring to (3). For the edge e with the ending vertex coordinates (x′,y′) and (x′′,y′′), its outward unit normal vector can be computed by nx=y′-y′′ and ny=x′′-x′.Assemble the coefficient matrices in (34) into the stiffness matrix of the global equilibrium equations (22) according to the rules below:(37a) K i i ⟵ - 1 2 ∫ e T i T d l N e C B - 1 2 B T C T N e T ∫ e T i T d l + δ p G e ∫ e T i T T i d l , (37b) K i j ⟵ - 1 2 ∫ e T i T d l N e C B + 1 2 B T C T N e T ∫ e T j T d l + δ p G e ∫ e T i T T j d l , (37c) K j i ⟵ - 1 2 ∫ e T j T d l N e C B + 1 2 B T C T N e T ∫ e T i T d l + δ p G e ∫ e T j T T i d l , (37d) K j j ⟵ 1 2 ∫ e T j T d l N e C B + 1 2 B T C T N e T ∫ e T j T d l + δ p G e ∫ e T j T T j d l ,where arrows indicate the operation to assemble submatrices into the stiffness matrix of the global equilibrium equations (22). The integrals in this section and Sections 3.5 and 3.6 can be exactly calculated by using the parametric equation of the edge e and the integral transform technique.


### 3.5. Displacement Boundary Constraint

Assume the exterior edgee of the element Ei is prescribed with the displacement boundary constraint u¯=u¯x,u¯yT. We can expand (15) as follows by using the notations of jump and average operators defined in (10) and (11):(38)as=-∫eσuine·vidl-∫eσvine·uidl+δpGe∫eui·vidl.Substituting (19) into (38) gives(39)as=-DiT∫eTiTdlNeCBDi-DiTBTCTNeT∫eTiTdlDi+δpGeDiT∫eTiTTidlDi,where Ne is the unit outward normal matrix introduced in (35). Assemble the coefficient matrices in (39) into the submatrix Kii in the stiffness matrix of (22):(40)-∫eTiTdlNeCB-BTCTNeT∫eTiTdl+δpGe∫eTiTTidl⟶Kii.Besides, substituting (19) into (18) leads to(41)lu=DiTBTCTNeT∫eu¯xu¯ydl+δpGeDiT∫eTiTu¯xu¯ydl,where the coefficient vectors are assembled to the right hand side Fi of the global equilibrium equations (22); that is,(42)BTCTNeT∫eu¯xu¯ydl+δpGe∫eTiTu¯xu¯ydl⟶Fi.


### 3.6. Surface Traction

Assume the elementEi subjected to the traction t¯=t¯x,t¯yT on its exterior edge e. Then, substituting (19) into (17) leads to(43)lt=DiT∫eTiTt¯xt¯ydl,where the integral on the right is assembled to the submatrix Fi in the right hand side of (22).


### 3.7. Contact Submatrices

In the IPG-DDA, blocks are subdivided and their boundaries become the combination of exterior element edges. Consequently, contacts need to be transferred to the element edges that compose original block boundaries. Contacts between elements may happen in three modes, that is, vertex-to-vertex, edge-to-edge, and vertex-to-edge cases, but the first two can always be decomposed into vertex-to-edge contacts. As defined in Section2.2, exterior element edges are portions of block boundaries, while interior elements edges are the others. We use the sign ∂Ei to indicate the set of the boundary of the element Ei.As shown in Figure3(a), assume that the vertex P1∈∂Ei contacts with the edge P2P3¯⊂∂Ej. A normal contact spring with the stiffness κ⊥ and a shear contact spring with the stiffness κ∥ are applied at the contact point P0∈∂Ej. Contact springs contribute the strain energies Π⊥ and Π∥:(44)Π⊥=12κ⊥h⊥2,(45)Π∥=12κ∥h∥2,where h⊥ and h∥ are the normal and the shear contact displacement, respectively, as illustrated in Figure 3(b). Based on the minimum potential energy principle, we need to minimize the strain energies Π⊥ and Π∥. To this end, h⊥ and h∥ have to be figured out in the first step.Figure 3
Contact between two elementsEi and Ej. (a) Implementation of normal contact spring and shear contact spring. (b) The normal contact distance h⊥ and the shear contact distance h∥.
(a) (b)Assume the coordinates of the pointPξ at the beginning and the end of the current time step are x¯k,y¯k and xk,yk, respectively, ξ=0,1,2,3. From Figure 3(b), we can find the fact that h⊥ is the distance between P1 and P2P3¯, and therefore we have(46)h⊥=Jd,where J denotes the area of the triangle of P1, P2, and P3 and d indicates the length of the edge P2P3¯:(47)J=1x1y11x2y21x3y3=1x¯1y¯11x¯2y¯21x¯3y¯3+1u1y¯11u2y¯21u3y¯3+1x¯1v11x¯2v21x¯3v3+1u1v11u2v21u3v3,(48)d=x2-x32+y2-y32.If we choose small enough time step length, the last item on the right side of (47) can be ignored, and d can also be approximated by the initial length of P2P3¯ during this time step:(49)J≈1x¯1y¯11x¯2y¯21x¯3y¯3+1u1y¯11u2y¯21u3y¯3+1x¯1v11y¯2v21y¯3v3≡J0+1u1y¯11u2y¯21u3y¯3+1x¯1v11y¯2v21y¯3v3,(50)d≈d0=x¯2-x¯32+y¯2-y¯32.Substitute (21), (49), and (50) into (46), and we have the following approximation of h⊥:(51)h⊥≈1d0J0+1d0y¯2-y¯3,x¯3-x¯2Tix1,y1Di+1d0y¯3-y¯1,x¯1-x¯3Tjx2,y2Dj+1d0y¯1-y¯2,x¯2-x¯1Tjx3,y3Dj.The shear contact displacementh∥ is actually the projection of the line P0P1¯ onto the edge P2P3¯, and therefore it can be determined as follows:(52)h∥=P0P1→·P2P3→P2P3¯=1dx3-x2,y3-y2x1-x0y1-y0.Due to the fact uξ=xξ-x¯ξ and vξ=yξ-y¯ξ, ξ=0,1,…,3, (52) becomes(53)h∥=1dx¯3-x¯2,y¯3-y¯2x¯1-x¯0y¯1-y¯0+1dx¯3-x¯2,y¯3-y¯2u1-u0v1-v0+1du3-u2,v3-v2x¯1-x¯0y¯1-y¯0+1du3-u2,v3-v2u1-u0v1-v0.When we use small enough time step length, the last item on the right side of the above equality can be ignored, and d can be approximated by d0. Then substituting (21) into (53) gives(54)h∥=B0d0-1+d0-1x¯3-x¯2,y¯3-y¯2Tix1,y1Di+d0-1x¯2-x¯3,y¯2-y¯3Tjx0,y0Dj,where B0=x¯1-x¯0x¯3-x¯2+y¯1-y¯0y¯3-y¯2.After substituting (51) into (44) and (54) into (45), we minimize the summation of Π⊥ and Π∥ with respect to Di and Dj and have the following contact submatrices:(55a)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Φ⊥TΦ⊥+κ∥Φ∥TΦ∥⟶Kii,(55b)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Φ⊥TΨ⊥+κ∥Φ∥TΨ∥⟶Kij,(55c)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Ψ⊥TΦ⊥+κ∥Ψ∥TΦ∥⟶Kji,(55d)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Ψ⊥TΨ⊥+κ∥Ψ∥TΨ∥⟶Kjj,(55e)-∂Π⊥0+Π∥0∂dri=κ⊥J0Φ⊥+κ∥B0Φ∥⟶Fi,(55f)-∂Π⊥0+Π∥0∂dsj=κ⊥J0Φ⊥+κ∥B0Φ∥⟶Fj,where r,s=1,2,…,6, dri and dsj denote the rth component of Di and the sth component of Dj, respectively, and the 1×6 matrices(56)Φ⊥=d0-1y¯2-y¯3,x¯3-x¯2Tix1,y1,Ψ⊥=d0-1y¯3-y¯1,x¯1-x¯3Tjx2,y2+d0-1y¯1-y¯2,x¯2-x¯1Tjx3,y3,Φ∥=d0-1y¯2-y¯3,x¯3-x¯2Tix1,y1,Ψ∥=d0-1x¯2-x¯3,y¯2-y¯3Tjx0,y0.The submatrices in (55a)–(55d) are assembled to the stiffness matrix of (22), while the submatrices in (55e) and (55f) are assembled to the right hand side of (22).


### 3.8. Friction Submatrices

Relative sliding occurs when the Mohr-Coulomb failure criterion is violated. In this case, normal contact springs need to be retained to prevent penetrations, but shear contact springs should be removed. For nonzero frictional angle, frictions have to be taken into account.As shown in Figure4(a), assume that the vertex P1∈∂Ei is sliding along the edge P2P3¯⊂Ej, and the normal contact spring with the stiffness κ⊥ is applied at the contact point P0∈P2P3¯. As shown in Figure 4(b), s0 and s1 are the sliding displacements of P0 and P1 along P2P3→, respectively, and h⊥ denotes the normal contact displacement. Then the normal contact force fn and the friction force f^ are(57)fn=κ⊥h⊥,(58)f^=sgn⁡P0P1→·P2P3→fntan⁡φ,where sgn⁡· is the sign function to determine the relative sliding direction, and φ>0 is the friction angle. The friction force f^ acts on both sides of the sliding surface. Consequently, the potential energy Πf contributed by the friction force f^ consists of the energy f^s0 for Ei and the energy f^s1 for Ej:(59)Πf=f^s0+f^s1.Based on the fact that s0 is the projection of the displacement of P0 on the contact surface P2P3→, while s1 is the projection of the movement of P1 on P2P3→, we can compute s0 and s1 as follows:(60)s0=-u0,v0·P2P3→P2P3¯=-u0,v0·1dx¯3-x¯2y¯3-y¯2,(61)s1=u1,v1·P2P3→P2P3¯=u1,v1·1dx¯3-x¯2y¯3-y¯2,where d is identical to the expression in (48) and can be approximated by d0 defined in (50) when using small enough time step length. Then, substitute (19) and (21) into (60) and (61), respectively, and we have(62)s0=-1d0DjTTjx0,y0Tx¯3-x¯2y¯3-y¯2,s1=1d0DiTTix1,y1Tx¯3-x¯2y¯3-y¯2.Now, using the expressions of s0 and s1 in the above two equalities, we minimize Πf with respect to Di and Dj, respectively:(63)-∂Πf0∂dri=-f^d0Tix1,y1Tx¯2-x¯3y¯2-y¯3⟶Fi,-∂Πf0∂dsj=f^d0Tjx0,y0Tx¯3-x¯2y¯3-y¯2⟶Fj,where r,s=1,2,…,6 and dri and dsj denote the rth component of Di and the sth component of Dj, respectively. The two vectors above are assembled to the right hand side of the global equilibrium equations (22).Figure 4
The vertexP1∈∂Ei slides along the edge P2P3¯⊂Ej. (a) Implementation of normal contact spring. (b) Sliding displacements s0 and s1.
(a) (b)


## 3.1. Inertia Force

In the time step intervalt,t+Δt, substituting (19) into (13) leads to(23)aρ=DiT∫EiρTiTTidE∂2∂t2Di≡DiTMi∂2∂t2Di,where the area integral in the parentheses can be calculated by the simplex integral formulations [39]. With the help of Newmark time integral scheme [40], we can solve the second derivative of Di in (23) as(24)aρ=1βΔt2DiTMiDi-1βΔtDiTMiVi0-12β-1DiTMiAi0,where γ and β indicate Newmark parameters and Vi0 and Ai0 denote the initial velocity and acceleration of the current time step. Notice that the initial velocity and acceleration in one time step are actually the finial velocity and acceleration in the last time step. Therefore we can use the following recursive formulations to update Vi0 and Ai0 in every time step:(25)Ai1=1βΔt2Di-1βΔtVi0-12β-1Ai0,Vi1=γβΔtDi+1-γβVi0+1-γ2βΔtAi0,where Vi1 and Ai1 refer to the finial velocity and acceleration, respectively. In the numerical examples of this article, we choose γ=1 and β=0.5 to avoid numerical damping and to ensure this time integration method to be unconditionally stable.Since the IPG method is based on the Galerkin weighted residual method, the global equilibrium equations have to hold for all test functionsvi, and therefore the vector DiT in (24) can be divided out from the equations. So we assemble 1/βΔt2Mi on the right side of (24) to [Kii] in the stiffness matrix of the global equilibrium equations (22), while -1/βΔtMiVi0 and -1/2β-1MiAi0 are added to Fi of the right hand side of (22). The submatrices derived below are assembled to the global equilibrium equations due to the same reason.


## 3.2. Element Stiffness

Denote the vector forms of the stress tensor and the strain tensor on the elementEi by σi and εi, respectively. Substituting (21) into (4) and (5), respectively, leads to(26)εi=εxεyγxy=000100000010000001Di≡BDi,σi=σxσyτxy=E¯1-ν¯21ν¯0ν¯10001-ν¯2εi≡CiBDi,where E¯=E and ν¯=ν for plane stress problems and E¯=E/(1-ν2) and ν¯=ν/(1-ν) for plane strain problems, in which the constants E and ν indicate Young’s modulus and Poisson’s ratio of the element Ei, respectively. The matrix [Ci] above is known as the material matrix of Ei. For the three-dimensional case rather than the two-dimensional one as concerned in this paper, contact points, edges, and planes have to be taken into account and detailed in different manners.Then substitute (19) and (26) into (14), and we have(27)aσ=∫EiDiTBTCiBDidEi=DiTEiBTCiBDi,where Ei is the area of the element Ei. The matrix EiBTCiB is assembled to the 6×6 submatrix Kii in the stiffness matrix of the global equilibrium equations (22).


## 3.3. Body Force

Assume the elementEi bearing the body force bi=bx,byT. Substituting (19) into (16) gives(28)lb=DiT∫EiTTbxbydEi,where the integral on the right of the above equality is assembled to the right hand side Fi of the global simultaneous equilibrium equations (22).


## 3.4. Interior Interface Stiffness

Assume the interior edgee as the common boundary of the elements Ei and Ej. Rewrite as in (15) into the summation below:(29)as=as,1+as,2+as,3,where(30)as,1=-∫eσune·vdl,as,2=-∫eσvne·udl,as,3=δpGe∫eu·vdl.Using the notations of average and jump operators defined in (10) and (11), respectively, we can expand the expressions above into the formulations below:(31)as,1=-12∫eσuine·vidl+12∫eσujne·vidl+12∫eσuine·vjdl-12∫eσujne·vjdl,(32)as,2=-12∫eσvine·uidl+12∫eσvjne·uidl+12∫eσvine·ujdl-12∫eσvjne·ujdl,(33)as,3=δpGe∫eui·vidl-δpGe∫eui·vjdl-δpGe∫euj·vidl+δpGe∫euj·vjdl.Then substitute (19) into (31), (32), and (33), respectively, and we have(34)as,1=-12DiT∫eTiTdlNeCBDi+12DjT∫eTiTdlNeCBDj+12DiT∫eTjTdlNeCBDi-12DjT∫eTjTdlNeCBDj,as,2=-12DiTBTCTNeT∫eTiTdlDi+12DjTBTCTNeT∫eTidlDi+12DiTBTCTNeT∫eTjTdlDj-12DjTBTCTNeT∫eTjTdlDj,as,3=δpGeDiT∫eTiTTidlDi-δpGeDjT∫eTjTTidlDi-δpGeDiT∫eTiTTjdlDj+δpGeDjT∫eTjTTjdlDj,where the matrix Ne is spanned by the unit outward normal vector nx,ny to the edge e:(35)Ne=nx0ny0nynx.Actually, Ne transfers the surface traction f¯x,f¯yT to the combination of stress components in the following way:(36)f¯xf¯y=σne=σxτxyτxyσynxny=Neσxσyτxy,where σ denotes the stress tensor referring to (3). For the edge e with the ending vertex coordinates (x′,y′) and (x′′,y′′), its outward unit normal vector can be computed by nx=y′-y′′ and ny=x′′-x′.Assemble the coefficient matrices in (34) into the stiffness matrix of the global equilibrium equations (22) according to the rules below:(37a) K i i ⟵ - 1 2 ∫ e T i T d l N e C B - 1 2 B T C T N e T ∫ e T i T d l + δ p G e ∫ e T i T T i d l , (37b) K i j ⟵ - 1 2 ∫ e T i T d l N e C B + 1 2 B T C T N e T ∫ e T j T d l + δ p G e ∫ e T i T T j d l , (37c) K j i ⟵ - 1 2 ∫ e T j T d l N e C B + 1 2 B T C T N e T ∫ e T i T d l + δ p G e ∫ e T j T T i d l , (37d) K j j ⟵ 1 2 ∫ e T j T d l N e C B + 1 2 B T C T N e T ∫ e T j T d l + δ p G e ∫ e T j T T j d l ,where arrows indicate the operation to assemble submatrices into the stiffness matrix of the global equilibrium equations (22). The integrals in this section and Sections 3.5 and 3.6 can be exactly calculated by using the parametric equation of the edge e and the integral transform technique.


## 3.5. Displacement Boundary Constraint

Assume the exterior edgee of the element Ei is prescribed with the displacement boundary constraint u¯=u¯x,u¯yT. We can expand (15) as follows by using the notations of jump and average operators defined in (10) and (11):(38)as=-∫eσuine·vidl-∫eσvine·uidl+δpGe∫eui·vidl.Substituting (19) into (38) gives(39)as=-DiT∫eTiTdlNeCBDi-DiTBTCTNeT∫eTiTdlDi+δpGeDiT∫eTiTTidlDi,where Ne is the unit outward normal matrix introduced in (35). Assemble the coefficient matrices in (39) into the submatrix Kii in the stiffness matrix of (22):(40)-∫eTiTdlNeCB-BTCTNeT∫eTiTdl+δpGe∫eTiTTidl⟶Kii.Besides, substituting (19) into (18) leads to(41)lu=DiTBTCTNeT∫eu¯xu¯ydl+δpGeDiT∫eTiTu¯xu¯ydl,where the coefficient vectors are assembled to the right hand side Fi of the global equilibrium equations (22); that is,(42)BTCTNeT∫eu¯xu¯ydl+δpGe∫eTiTu¯xu¯ydl⟶Fi.


## 3.6. Surface Traction

Assume the elementEi subjected to the traction t¯=t¯x,t¯yT on its exterior edge e. Then, substituting (19) into (17) leads to(43)lt=DiT∫eTiTt¯xt¯ydl,where the integral on the right is assembled to the submatrix Fi in the right hand side of (22).


## 3.7. Contact Submatrices

In the IPG-DDA, blocks are subdivided and their boundaries become the combination of exterior element edges. Consequently, contacts need to be transferred to the element edges that compose original block boundaries. Contacts between elements may happen in three modes, that is, vertex-to-vertex, edge-to-edge, and vertex-to-edge cases, but the first two can always be decomposed into vertex-to-edge contacts. As defined in Section2.2, exterior element edges are portions of block boundaries, while interior elements edges are the others. We use the sign ∂Ei to indicate the set of the boundary of the element Ei.As shown in Figure3(a), assume that the vertex P1∈∂Ei contacts with the edge P2P3¯⊂∂Ej. A normal contact spring with the stiffness κ⊥ and a shear contact spring with the stiffness κ∥ are applied at the contact point P0∈∂Ej. Contact springs contribute the strain energies Π⊥ and Π∥:(44)Π⊥=12κ⊥h⊥2,(45)Π∥=12κ∥h∥2,where h⊥ and h∥ are the normal and the shear contact displacement, respectively, as illustrated in Figure 3(b). Based on the minimum potential energy principle, we need to minimize the strain energies Π⊥ and Π∥. To this end, h⊥ and h∥ have to be figured out in the first step.Figure 3
Contact between two elementsEi and Ej. (a) Implementation of normal contact spring and shear contact spring. (b) The normal contact distance h⊥ and the shear contact distance h∥.
(a) (b)Assume the coordinates of the pointPξ at the beginning and the end of the current time step are x¯k,y¯k and xk,yk, respectively, ξ=0,1,2,3. From Figure 3(b), we can find the fact that h⊥ is the distance between P1 and P2P3¯, and therefore we have(46)h⊥=Jd,where J denotes the area of the triangle of P1, P2, and P3 and d indicates the length of the edge P2P3¯:(47)J=1x1y11x2y21x3y3=1x¯1y¯11x¯2y¯21x¯3y¯3+1u1y¯11u2y¯21u3y¯3+1x¯1v11x¯2v21x¯3v3+1u1v11u2v21u3v3,(48)d=x2-x32+y2-y32.If we choose small enough time step length, the last item on the right side of (47) can be ignored, and d can also be approximated by the initial length of P2P3¯ during this time step:(49)J≈1x¯1y¯11x¯2y¯21x¯3y¯3+1u1y¯11u2y¯21u3y¯3+1x¯1v11y¯2v21y¯3v3≡J0+1u1y¯11u2y¯21u3y¯3+1x¯1v11y¯2v21y¯3v3,(50)d≈d0=x¯2-x¯32+y¯2-y¯32.Substitute (21), (49), and (50) into (46), and we have the following approximation of h⊥:(51)h⊥≈1d0J0+1d0y¯2-y¯3,x¯3-x¯2Tix1,y1Di+1d0y¯3-y¯1,x¯1-x¯3Tjx2,y2Dj+1d0y¯1-y¯2,x¯2-x¯1Tjx3,y3Dj.The shear contact displacementh∥ is actually the projection of the line P0P1¯ onto the edge P2P3¯, and therefore it can be determined as follows:(52)h∥=P0P1→·P2P3→P2P3¯=1dx3-x2,y3-y2x1-x0y1-y0.Due to the fact uξ=xξ-x¯ξ and vξ=yξ-y¯ξ, ξ=0,1,…,3, (52) becomes(53)h∥=1dx¯3-x¯2,y¯3-y¯2x¯1-x¯0y¯1-y¯0+1dx¯3-x¯2,y¯3-y¯2u1-u0v1-v0+1du3-u2,v3-v2x¯1-x¯0y¯1-y¯0+1du3-u2,v3-v2u1-u0v1-v0.When we use small enough time step length, the last item on the right side of the above equality can be ignored, and d can be approximated by d0. Then substituting (21) into (53) gives(54)h∥=B0d0-1+d0-1x¯3-x¯2,y¯3-y¯2Tix1,y1Di+d0-1x¯2-x¯3,y¯2-y¯3Tjx0,y0Dj,where B0=x¯1-x¯0x¯3-x¯2+y¯1-y¯0y¯3-y¯2.After substituting (51) into (44) and (54) into (45), we minimize the summation of Π⊥ and Π∥ with respect to Di and Dj and have the following contact submatrices:(55a)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Φ⊥TΦ⊥+κ∥Φ∥TΦ∥⟶Kii,(55b)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Φ⊥TΨ⊥+κ∥Φ∥TΨ∥⟶Kij,(55c)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Ψ⊥TΦ⊥+κ∥Ψ∥TΦ∥⟶Kji,(55d)∂2Π⊥+Π∥∂dri∂dsj=κ⊥Ψ⊥TΨ⊥+κ∥Ψ∥TΨ∥⟶Kjj,(55e)-∂Π⊥0+Π∥0∂dri=κ⊥J0Φ⊥+κ∥B0Φ∥⟶Fi,(55f)-∂Π⊥0+Π∥0∂dsj=κ⊥J0Φ⊥+κ∥B0Φ∥⟶Fj,where r,s=1,2,…,6, dri and dsj denote the rth component of Di and the sth component of Dj, respectively, and the 1×6 matrices(56)Φ⊥=d0-1y¯2-y¯3,x¯3-x¯2Tix1,y1,Ψ⊥=d0-1y¯3-y¯1,x¯1-x¯3Tjx2,y2+d0-1y¯1-y¯2,x¯2-x¯1Tjx3,y3,Φ∥=d0-1y¯2-y¯3,x¯3-x¯2Tix1,y1,Ψ∥=d0-1x¯2-x¯3,y¯2-y¯3Tjx0,y0.The submatrices in (55a)–(55d) are assembled to the stiffness matrix of (22), while the submatrices in (55e) and (55f) are assembled to the right hand side of (22).


## 3.8. Friction Submatrices

Relative sliding occurs when the Mohr-Coulomb failure criterion is violated. In this case, normal contact springs need to be retained to prevent penetrations, but shear contact springs should be removed. For nonzero frictional angle, frictions have to be taken into account.As shown in Figure4(a), assume that the vertex P1∈∂Ei is sliding along the edge P2P3¯⊂Ej, and the normal contact spring with the stiffness κ⊥ is applied at the contact point P0∈P2P3¯. As shown in Figure 4(b), s0 and s1 are the sliding displacements of P0 and P1 along P2P3→, respectively, and h⊥ denotes the normal contact displacement. Then the normal contact force fn and the friction force f^ are(57)fn=κ⊥h⊥,(58)f^=sgn⁡P0P1→·P2P3→fntan⁡φ,where sgn⁡· is the sign function to determine the relative sliding direction, and φ>0 is the friction angle. The friction force f^ acts on both sides of the sliding surface. Consequently, the potential energy Πf contributed by the friction force f^ consists of the energy f^s0 for Ei and the energy f^s1 for Ej:(59)Πf=f^s0+f^s1.Based on the fact that s0 is the projection of the displacement of P0 on the contact surface P2P3→, while s1 is the projection of the movement of P1 on P2P3→, we can compute s0 and s1 as follows:(60)s0=-u0,v0·P2P3→P2P3¯=-u0,v0·1dx¯3-x¯2y¯3-y¯2,(61)s1=u1,v1·P2P3→P2P3¯=u1,v1·1dx¯3-x¯2y¯3-y¯2,where d is identical to the expression in (48) and can be approximated by d0 defined in (50) when using small enough time step length. Then, substitute (19) and (21) into (60) and (61), respectively, and we have(62)s0=-1d0DjTTjx0,y0Tx¯3-x¯2y¯3-y¯2,s1=1d0DiTTix1,y1Tx¯3-x¯2y¯3-y¯2.Now, using the expressions of s0 and s1 in the above two equalities, we minimize Πf with respect to Di and Dj, respectively:(63)-∂Πf0∂dri=-f^d0Tix1,y1Tx¯2-x¯3y¯2-y¯3⟶Fi,-∂Πf0∂dsj=f^d0Tjx0,y0Tx¯3-x¯2y¯3-y¯2⟶Fj,where r,s=1,2,…,6 and dri and dsj denote the rth component of Di and the sth component of Dj, respectively. The two vectors above are assembled to the right hand side of the global equilibrium equations (22).Figure 4
The vertexP1∈∂Ei slides along the edge P2P3¯⊂Ej. (a) Implementation of normal contact spring. (b) Sliding displacements s0 and s1.
(a) (b)


## 4. Implementation Aspects

### 4.1. Solution Algorithm

The computer program of the proposed IPG-DDA has been coded based on the formulations derived in the previous section. The process of the program can be summarized as the following steps that consist of the time step loop2⇌5 and the open-close iteration 4⇌5:(1)
Initialize the time step countern←0; then input data and set the maximum time step number n∗.(2)
n←n+1,stop the program if n>n∗, and otherwise detect contacts and determine their types.(3)
Assemble the simultaneous global equilibrium equations of the IPG-DDA.(4)
Solve the equations by block Jacobi preconditioned conjugate gradient (BJ-PCG) method [41].(5)
Figure out vertex displacements and contact distances and judge whether all contact springs satisfy no tension in an open contact and no penetration in a close contact:if true, then update element displacements, stresses, and velocities, andgo to (2);if false, modify the contact springs which violate the criterion, and reassemble the global equilibrium equations; thengo to (4).


### 4.2. Auxiliary Criteria for Contact Detection

During an open-close iteration, the status and parameters of closed contacts must be saved and passed on to the next step. These pieces of information are quite different for vertex-to-vertex contacts and vertex-to-edge contacts, where every vertex-to-vertex contact with parallel edges is represented by two vertex-to-edge contacts. Therefore, it is important to determine the types of all contacts. The DDA carries out this task based on two criteria: the distances between vertexes and the degrees of overlapping angles. Nevertheless, it is inadequate to use them to determine the contacts in the IPG-DDA.Block boundaries become the broken lines composed of exterior element edges in our scheme, as shown in Figure1. Since contacts can only occur along block surfaces in reality, contact springs should not appear on interior element edges. Such redundant springs will lead to fictitious stress concentration phenomena and increase the loops of open-close iterations. For this reason, we introduce the following two auxiliary contact detection criteria:(C1)
The elements having no edge exposed to a contact interface do not participate the contact.(C2)
Interior element edges cannot be counted in as entrance lines.The first filters out the interior elements whose edges are all interior edges before using DDA contact detection criteria, while the second constraint deletes the contact pairs containing interior edges after employing all DDA contact detection criteria.Take the block contact in Figure1 as an example, where two blocks are in the surface-to-surface contact between the broken lines P1P2-P2+P3-P3+P4-P4+P5¯ and P6P7-P7+P8-P8+P9-P9+P10¯. The letter P is omitted in the figure. Here, Pi- and Pi+ are located on the same node but belong to different elements, i=2,3,…,9. Based on DDA contact detection criteria and (C1), elements contact between the edge couples: P1P2-¯ and P6P7-¯, P2+P3-¯ and P7+P8-¯, P3+P4-¯ and P8+P9-¯, and P4+P5¯ and P9+P10¯. These contacts consist of the vertex-to-edge contacts as enumerated in Table 1. Note that the vertex couples P2+ and P7-, P2- and P7+, P3+ and P8-, P3- and P8+, P4+ and P9-, and P4- and P9+ are not vertex-to-vertex contacts. Take P2+ and P7- as an example to explain the reasons. According to (C2), P2+ is not allowed contact with the interior edge P7-P11, while the situation that P2+ contacts with P7-P8+¯ has been taken into account as the contact of the vertex P2- with the edge P7-P8+¯.Table 1
Invasion vertexes and entrance lines of the contact pairs between Block 1 and Block 2 in Figure1.
Pair number Invasion vertex Entrance line Invasion vertex Entrance line 1 P 1 P 6 P 7 + ¯ P 6 P 1 P 2 + ¯ 2 P 2 + P 7 + P 6 ¯ P 7 + P 2 + P 1 ¯ 3 P 2 - P 7 - P 8 + ¯ P 7 - P 2 - P 3 + ¯ 4 P 3 + P 8 + P 7 - ¯ P 8 + P 3 + P 2 - ¯ 5 P 3 - P 8 - P 9 + ¯ P 8 - P 3 - P 4 + ¯ 6 P 4 + P 9 + P 8 - ¯ P 9 + P 4 + P 3 - ¯ 7 P 4 - P 9 - P 10 ¯ P 9 - P 4 - P 5 ¯ 8 P 5 P 10 P 9 - ¯ P 10 P 5 P 4 - ¯This example demonstrates the merit of the auxiliary criteria (C1) and (C2) that a lot of computational cost can be saved in the phase for determining contact types. When using very fine meshes on contact zones, this merit will be much more salient to guarantee convergent open-close iteration and to improve the computing efficiency.


## 4.1. Solution Algorithm

The computer program of the proposed IPG-DDA has been coded based on the formulations derived in the previous section. The process of the program can be summarized as the following steps that consist of the time step loop2⇌5 and the open-close iteration 4⇌5:(1)
Initialize the time step countern←0; then input data and set the maximum time step number n∗.(2)
n←n+1,stop the program if n>n∗, and otherwise detect contacts and determine their types.(3)
Assemble the simultaneous global equilibrium equations of the IPG-DDA.(4)
Solve the equations by block Jacobi preconditioned conjugate gradient (BJ-PCG) method [41].(5)
Figure out vertex displacements and contact distances and judge whether all contact springs satisfy no tension in an open contact and no penetration in a close contact:if true, then update element displacements, stresses, and velocities, andgo to (2);if false, modify the contact springs which violate the criterion, and reassemble the global equilibrium equations; thengo to (4).


## 4.2. Auxiliary Criteria for Contact Detection

During an open-close iteration, the status and parameters of closed contacts must be saved and passed on to the next step. These pieces of information are quite different for vertex-to-vertex contacts and vertex-to-edge contacts, where every vertex-to-vertex contact with parallel edges is represented by two vertex-to-edge contacts. Therefore, it is important to determine the types of all contacts. The DDA carries out this task based on two criteria: the distances between vertexes and the degrees of overlapping angles. Nevertheless, it is inadequate to use them to determine the contacts in the IPG-DDA.Block boundaries become the broken lines composed of exterior element edges in our scheme, as shown in Figure1. Since contacts can only occur along block surfaces in reality, contact springs should not appear on interior element edges. Such redundant springs will lead to fictitious stress concentration phenomena and increase the loops of open-close iterations. For this reason, we introduce the following two auxiliary contact detection criteria:(C1)
The elements having no edge exposed to a contact interface do not participate the contact.(C2)
Interior element edges cannot be counted in as entrance lines.The first filters out the interior elements whose edges are all interior edges before using DDA contact detection criteria, while the second constraint deletes the contact pairs containing interior edges after employing all DDA contact detection criteria.Take the block contact in Figure1 as an example, where two blocks are in the surface-to-surface contact between the broken lines P1P2-P2+P3-P3+P4-P4+P5¯ and P6P7-P7+P8-P8+P9-P9+P10¯. The letter P is omitted in the figure. Here, Pi- and Pi+ are located on the same node but belong to different elements, i=2,3,…,9. Based on DDA contact detection criteria and (C1), elements contact between the edge couples: P1P2-¯ and P6P7-¯, P2+P3-¯ and P7+P8-¯, P3+P4-¯ and P8+P9-¯, and P4+P5¯ and P9+P10¯. These contacts consist of the vertex-to-edge contacts as enumerated in Table 1. Note that the vertex couples P2+ and P7-, P2- and P7+, P3+ and P8-, P3- and P8+, P4+ and P9-, and P4- and P9+ are not vertex-to-vertex contacts. Take P2+ and P7- as an example to explain the reasons. According to (C2), P2+ is not allowed contact with the interior edge P7-P11, while the situation that P2+ contacts with P7-P8+¯ has been taken into account as the contact of the vertex P2- with the edge P7-P8+¯.Table 1
Invasion vertexes and entrance lines of the contact pairs between Block 1 and Block 2 in Figure1.
Pair number Invasion vertex Entrance line Invasion vertex Entrance line 1 P 1 P 6 P 7 + ¯ P 6 P 1 P 2 + ¯ 2 P 2 + P 7 + P 6 ¯ P 7 + P 2 + P 1 ¯ 3 P 2 - P 7 - P 8 + ¯ P 7 - P 2 - P 3 + ¯ 4 P 3 + P 8 + P 7 - ¯ P 8 + P 3 + P 2 - ¯ 5 P 3 - P 8 - P 9 + ¯ P 8 - P 3 - P 4 + ¯ 6 P 4 + P 9 + P 8 - ¯ P 9 + P 4 + P 3 - ¯ 7 P 4 - P 9 - P 10 ¯ P 9 - P 4 - P 5 ¯ 8 P 5 P 10 P 9 - ¯ P 10 P 5 P 4 - ¯This example demonstrates the merit of the auxiliary criteria (C1) and (C2) that a lot of computational cost can be saved in the phase for determining contact types. When using very fine meshes on contact zones, this merit will be much more salient to guarantee convergent open-close iteration and to improve the computing efficiency.


## 5. Numerical Examples

In this section, three numerical examples on contact problems are conducted and solved by the IPG-DDA. The results are also compared with analytical solutions and ANSYS simulations.Example 1.
The first numerical example is designed to examine the validity of the IPG-DDA for contact interactions. Consider the problem as illustrated in Figure5(a), where the linearly elastic Block 1 rests on the rigid Block 2. The top surface of Block 1 is subjected to a uniform pressure with the intensity of 1 MPa. The displacements of the bottom and the two sides of Block 2 are completely constrained. The geometry sizes and the coordinate system are annotated in Figure 5(a), and the material properties of the two blocks are given in Table 2. According to the loading condition, we can easily know the analytical contact pressure is −1 MPa distributed along the contact region AB¯ as shown in Figure 5(a).
We apply both ANSYS and the IPG-DDA to solve this example problem. ANSYS gives reference results. In both of their computational models, only the right half of the problem is taken into account owing to the symmetry. Block 1 is discretized by the mesh as shown in Figure5(b). But the treatments for Block 2 are different: it is simulated by a pilot node in ANSYS (a pilot node is an element with only one node, which is provided to simulate rigid target surfaces in ANSYS), while it is modeled as an integrated block without subdivision in the IPG-DDA. Detailed computational settings are stated in Table 3. It can be noticed that the total DOFs required by the IPG-DDA are only half that by ANSYS with the same mesh. During the computation of the IPG-DDA, we adopt δp=5 for the penalties appearing in (37a)–(37d), (40), and (42).
Comparative Task 1. Firstly, we compare the results from the IPG-DDA and ANSYS on the contact pressures distributed along the right half of the interface AB¯. As shown in Figure 6, it can be observed that both of the solutions computed by the IPG-DDA and ANSYS coincide with the analytical contact pressure curve. The absolute errors of the two methods are all less than 1e10.
Comparative Task 2. The results acquired from the IPG-DDA and ANSYS are compared on the horizontal displacement component ux of the right half of Block 1. Figures 7(a) and 7(b) draw the contour plots of ux computed by the two methods, respectively. It is shown that the displacement varying in the x direction of Block 1 figured out by our scheme is in good accordance with the simulation of ANSYS.
Comparative Task 3. Then we compare the solutions of the IPG-DDA and ANSYS for the vertical displacement component uy of the right half of Block 1. The contour plots of uy computed by the two methods are illustrated in Figures 8(a) and 8(b), respectively. From the figures, we can find that the result of the IPG-DDA agrees very well with that of ANSYS for the displacement uy of Block 1.Table 2
Material property parameters for Block 1 and Block 2 in Example1.
⁢ Block 1 Block 2 Plane stress Yes Yes Mass density (Kg/m2) 0 0 Young’s modulus (GPa) 10 1000 Poisson’s ratio 0.3 0Table 3
Computational settings used in ANSYS and the IPG-DDA for Example1.
⁢ ANSYS IPG-DDA Element for Block 1 PLANE82 Constant stress element Element for Block 2 Pilot node DDA block Contact element CONTA172 and TARGE169 — DOFs per element 12/PLANE82 and 3/pilot 6 Number of elements 888 PLANE82s and 1 pilot 888 elements and 1 block Total DOFs 10659 5334 Normal contact stiffness (N/m) 1E14 1E14 Shear contact stiffness (N/m) 3 E 13 † 3E13 †In Examples 1, 2, and 3, the shear contact stiffness takes the magnitude of Poisson’s ratio times the normal contact stiffness.Figure 5
The linearly elastic Block 1 rests on the rigid Block 2.
(a)
Geometry sizes and the coordinate system (b)
Mesh for the right half of Block 1Figure 6
The computed contact pressures along the right half of the interfaceAB¯.Figure 7
Displacements (m) in thex-direction obtained by different methods.
(a)
IPG-DDA result (b)
ANSYS resultFigure 8
Displacements (m) in they-direction obtained by different methods.
(a)
IPG-DDA result (b)
ANSYS resultExample 2.
This numerical example aims to verify the accuracy and effectiveness of the IPG-DDA for flexible-flexible contacts. Consider the Hertz contact problem as shown in Figure9(a), where the cylinder Block 1 is against the top surface of Block 2. Both of them are assumed to be perfectly elastic. And their material parameters are given in Table 4. Block 1 is subjected to the pressure P with the intensity of 1000 N/mm along its top generatrix; the displacement of the bottom of Block 2 is constrained. The geometry sizes and the coordinate system are illustrated in Figure 9(a). According to Hertz contact theory [42], the half width b of the contact region is(64)b=1.522PRE,where R is the radium of the cylinder and E is Young’s modulus; the normal contact pressure qx varying along the interface AB¯ has the following analytical solution:(65)qx=2Pπb2b2-x2,0⩽x⩽b,and the maximum contact pressure q0 takes place at the origin of the coordinate system:(66)q0=0.418PER.For Example 2, the theoretical half width of the contact region is b≈3.321276mm, and the analytical maximum contact pressure q0≈383.103328MPa.
We implement the IPG-DDA to solve this problem and use ANSYS to provide the reference results. Considering the symmetry existing in this example, we need to analyze only the right-lower quarter of Block 1 and the right half of Block 2. The mesh employed by the IPG-DDA and ANSYS is given in Figures9(b) and 9(c). The computational settings of the two methods are listed in Table 5. Besides, we adopt δp=100 for the penalties used in (37a)–(37d), (40), and (42) during the computation of the IPG-DDA.
Comparative Task 1. We first compare the contact pressures computed by the IPG-DDA and ANSYS with the theoretical solution. As shown in Figure 10, the normal contact pressure curves obtained from different methods are plotted. It can be found that both the simulations by the IPG-DDA and ANSYS are in good agreement with the analytical solution. From the figure, we can also observe that both of the maximum contact pressures obtained by the two methods appear at the contact center A with the coordinates (0,500). The IPG-DDA gives the peak normal contact pressure −383.508727 MPa, while ANSYS gives the value −383.485435 MPa. Compared with the analytical maximum contact pressure q0, the relative errors of the IPG-DDA and ANSYS are 0.11% and 0.10%, respectively. This illustrates that the proposed method is as accurate as ANSYS on the solution of the contact pressure distribution.
Comparative Task 2. Next, we compare the computed displacements of Block 2 restricted to the contact domain. Figure 11(a) plots the horizontal displacement component ux obtained by the IPG-DDA and ANSYS, while Figure 11(b) draws the vertical displacement component uy of the results from the two methods. It is shown in the two figures that the solution of the proposed method conforms very well with that of ANSYS.
Comparative Task 3. Then, we compare the displacement variations over the entire computing model solved by the IPG-DDA and ANSYS. Figures 12(a) and 12(b) provide the contour plots of the displacement distribution in the x direction from the two methods, respectively. Figures 12(c) and 12(d) present the contour images of the displacement in the y direction computed by the proposed coupling scheme and ANSYS, respectively. It can be observed that the results obtained by our method agree with the solutions given by ANSYS very well.
Comparative Task 4. Now, the variations of the stress component σx of both Block 1 and Block 2 solved by the IPG-DDA and ANSYS are compared. Figure 13(a) plots the contour distribution of the σx figured out by the proposed method, and Figure 13(b) zooms in to illustrate the details on the variable around the contact domain. Figure 13(c) supplies the contour picture of σx computed by ANSYS, and Figure 13(d) also displays the situation nearby the contact region. It is shown that the result of our method is in fine accordance with that of ANSYS.
Comparative Task 5. We bring the distributions of the stress component σy solved out by the IPG-DDA and ANSYS together to compare them. Figure 14(a) draws the contour of σy obtained by our coupling scheme, and Figure 14(b) enlarges the picture to show the detailed variation within the neighborhood of the contact. Figure 14(c) plots the contour of the σy given by ANSYS, and Figure 14(d) zooms in to display the stress around the contact domain. It can be found that the solution on σy of our method agrees with that of ANSYS very well. If it is noticed that the DOFs consumed by the IPG-DDA are only half of those by ANSYS, we can realize the high efficiency of the proposed method.Table 4
Material properties of Block 1 and Block 2 in Example2.
Block 1 Block 2 Plane stress Yes Yes Mass density 0 0 Young’s modulus 2.1E5 MPa 2.1E5 MPa Poisson’s ratio 0.3 0.3Table 5
Computational settings used in ANSYS and the IPG-DDA for Example2.
ANSYS IPG-DDA Element PLANE82 Constant stress element Contact element CONTA172 with TARGE169 — DOFs per element 12 6 Number of elements 7812 7812 Total DOFs 93744 46872 Normal contact stiffness (N/m) A L M † 2.3 E 8 ‡ Shear contact stiffness (N/m) ALM 6.9E7 Number of substeps 3 1 †The normal contact stiffness is determined by the augmented Lagrangian method in ANSYS. ‡The normal contact stiffness equals the equivalent Young’s modulus under the plane strain condition in the IPG-DDA.Figure 9
An elastic cylinder subjected to a uniform pressureP pressed against the surface of another elastic block.
(a)
Geometry sizes (b)
Triangular mesh (c)
Mesh detail nearby the contact regionFigure 10
Normal contact pressure (MPa) within the contact region computed by the IPG-DDA, ANSYS, and the analytical solution.Figure 11
The displacement of Block 2 changing within the contact region obtained by different methods.
(a)
Horizontal displacement componentux (b)
Vertical displacement componentuyFigure 12
The contours of the displacement componentsux (mm) and uy (mm) of the model obtained by different methods.
(a)
u x computed by the IPG-DDA (b)
u x computed by ANSYS (c)
u y computed by the IPG-DDA (d)
u y computed by ANSYSFigure 13
The contours of the stressσx (MPa) solved by different methods.
(a)
σ x computed by the IPG-DDA (b)
σ x near the contact domain simulated by the IPG-DDA (c)
σ x computed by ANSYS (d)
σ x near the contact domain simulated by ANSYSFigure 14
The contours of the stressσy (MPa) solved by different methods.
(a)
σ y computed by the IPG-DDA (b)
σ y near the contact domain simulated by the IPG-DDA (c)
σ y computed by ANSYS (d)
σ y near the contact domain simulated by ANSYSExample 3.
The third numerical example is to validate the capability of the IPG-DDA for simulating the dynamic contact problem of two deformable blocks in a relative sliding. As shown in Figure15(a), consider two linearly elastic trapezoidal blocks, Block 1 and Block 2. The computing time t is limited in the range [0,0.02] s. During the whole analysis, Block 1 is always subjected to a uniform pressure with the intensity q=1000 N/m distributed on its top surface. At the initial time t=0 s, Block 1 rests on the top of Block 2. The geometry sizes and the adopted coordinate system are marked in Figure 15(a), and the material properties of the two blocks are supplied in Table 6. In order to study the pure responses of Block 1 to external tractions, we do not take gravity acceleration into account.
We implement both the IPG-DDA and ANSYS to simulate the sliding process. Figure15(b) shows the mesh adopted in the two methods, and the parameters used in their computations are listed in Table 7. Besides, we adopt the penalty δp=1000 in (37a)–(37d), (40), and (42) for the IPG-DDA.
Comparative Task 1. Firstly, we compare the displacements of the point D varying in the time interval [0,0.02] obtained by the IPG-DDA and ANSYS. The point D located at coordinates (4,0.25) is one of the corners of Block 1. Figures 16(a) and 16(b) plot the computed displacement of D in the x and y directions, respectively. In the two figures, the analytical solutions are computed based on Newton’s laws and the theory of kinematics as follows: (67)ux=12·qLρS·cos⁡θ·sin⁡θ·t2,uy=12·qLρS·sin2θ·t2,0⩽t⩽0.02,where q=1000 N/m is the traction intensity, L=4m is the top surface length of Block 1, ρ=0.01Kg/m2 is the mass density, S=9 m2 is the area of Block 1, θ=arctan⁡(1/8) is the incline angle of the top surface of Block 2, and t indicates the time variable. As the time step lengths and the total time step numbers of the IPG-DDA and ANSYS are not identical, their results on the time nodes of themselves are not overlapped with each other in Figure 16. However, we can observe from the plot that both the results of our method and ANSYS are in good accordance with the analytical solution (67).
Comparative Task 2. Next, we compare the distributions of the stress component σx obtained by the IPG-DDA and ANSYS. Figures 17(a) and 17(b) show the contour plots of σx at t=0.01 s figured out by the two methods, respectively. Figures 17(c) and 17(d) illustrate the distributions of σx at t=0.02 s computed by our method and ANSYS, respectively. Studying these figures, we can find that the simulation results of the proposed coupling scheme agree with that of ANSYS very well.
Comparative Task 3. Then, the stress components σy solved out by the IPG-DDA and ANSYS are compared. Figures 18(a) and 18(b) give the contour distributions of σy at t=0.01 s obtained by the two methods, respectively. Figures 18(c) and 18(d) draw the contours of σy at t=0.02 s computed by the proposed scheme and ANSYS, respectively. It can be observed from the plots that the results of our method conform very well with that of ANSYS.
Comparative Task 4. Now, we take Coulomb’s friction into account for this example. As shown in Figure 15(a), it can be easily figured out that the interface between Block 1 and Block 2 is inclined at the angle θ with the tangent of tan⁡θ=0.125 to the negative direction of the x-axis. According to kinetics theory, a stick response is expected when the friction coefficient tan⁡φ used in (58) is larger than tan⁡θ. We run a couple of computations by using the IPG-DDA as follows: the first adopts the joint property with the friction coefficient tan⁡φ=0.12<tan⁡θ such that a slipping response is predicted along the interface; and the second chooses the joint property with tan⁡φ=0.126>tan⁡θ so that a stick state is expected.
Figure19 compares the computed geometry configurations of the considered blocks with the friction coefficients tan⁡φ=0.12 and 0.126 by the IPG-DDA. Figures 19(a) and 19(b) plot the case with the friction coefficient tan⁡φ=0.12 at the time instants t=0.01 s and t=0.02 s, respectively, while Figures 19(c) and 19(d) illustrate the case with the friction coefficient tan⁡φ=0.12 at t=0.01 s and t=0.02 s, respectively. In the first case, we have tan⁡φ<tan⁡θ, so a slipping response is expected. In the second case, there is tan⁡φ>tan⁡θ; thus the stick states are predicted. It can be observed from Figure 19 that these friction simulations by the IPG-DDA agree with the predictions based on kinetics very well.
In Figure20, the distributions of the block horizontal displacement components ux solved out by the IPG-DDA are compared in the cases of using different friction coefficients. Figures 20(a) and 20(b) show the contours of ux when tan⁡φ=0.12 at t=0.01 s and t=0.02 s, respectively. In this case, the slipping response is expected owing to tan⁡φ<tan⁡θ. The rigid relative movements between Block 1 and Block 2 along the inclined plane can be found in the two figures, while Figures 20(c) and 20(d) illustrate the contours of ux when tan⁡φ=0.126 at t=0.01 s and t=0.02 s, respectively. In this case, the stick state is predicted because of tan⁡φ>tan⁡θ. The continuously changing displacement fields of both Block 1 or Block 2 in the two figures reflect the distinguished effects due to the different sizes of the adopted friction coefficients. These observations from Figure 20 conform to the predictions.
Figure21 compares the distributions of the block vertical displacement components uy obtained by the IPG-DDA for different friction coefficients. Figures 21(a) and 21(b) plot the contours of uy for tan⁡φ=0.12<tan⁡θ at t=0.01 s and t=0.02 s, respectively, and a slipping response is expected in this situation of the fiction coefficient. Figures 21(c) and 21(d) draw the contours of uy for tan⁡φ=0.126>tan⁡θ at t=0.01 s and t=0.02 s, respectively, and the stick state is expected under this fictional condition. Distinctly discontinuous distribution of uy can be found along the interface of Block 1 and Block 2 as shown in Figures 21(a) and 21(b). But the distribution of uy can be observed to be nearly continuous along the block interface plane from Figures 21(c) and 21(d). Such different results between the two cases are in accordance with the predictions for the using of distinguished friction coefficients.Table 6
Material properties of Block 1 and Block 2 in Example3.
Block 1 Block 2 Plane strain Yes Yes Mass densityρ (Kg/m2) 0.01 0.01 Young’s modulus (Pa) 2.1E8 2.1E8 Poisson’s ratio 0.3 0.3 Gravity acceleration (m/s2) 0 0Table 7
Computational settings used in ANSYS and IPG-DDA for Example3.
ANSYS IPG-DDA Element type PLANE82 element Constant stress element Contact element type CONTA172 and TARGE169 — DOFs per element 12 6 Number of elements 2356 2356 Total DOFs 28272 14136 Normal contact stiffness (N/m) 1E9 1E9 Shear contact stiffness (N/m) 3E8 3E8 Time step length range (s) 5 E - 6,5 E - 5 † 5 E - 5,5 E - 4 ‡ †ANSYS command DELTIM is used to specify the time step sizes which are automatically adjusted to vary in the range. ‡The IPG-DDA inherits the open-close iteration algorithm from the DDA which reduces time step sizes once an open-close iteration process cannot reach the convergence criterion within six cycles.Figure 15
The dynamic relative sliding contact between two linearly elastic blocks along the plane with the inclination angleθ satisfying tan⁡θ=0.125 to the negative direction of the x-axis.
(a)
Problem geometry (b)
Mesh configurationFigure 16
The displacement of the pointD(4,2.5) obtained by different methods on sampled time nodes. The sampling time nodes are every 18 nodes in ANSYS results, while they are every 4 nodes in the IPG-DDA results.
(a)
The displacement componentux (m) along the ox-axis (b)
The displacement componentuy (m) along the oy-axisFigure 17
The stress componentσx (Pa) computed by different methods at t=0.01 s and t=0.02 s.
(a)
σ x computed by the IPG-DDA at t=0.01 (b)
σ x computed by ANSYS at t=0.01 (c)
σ x computed by the IPG-DDA at t=0.02 (d)
σ x computed by ANSYS at t=0.02Figure 18
The stress componentσy (Pa) computed by different methods at t=0.01 s and t=0.02 s.
(a)
σ y computed by the IPG-DDA at t=0.01 (b)
σ y computed by ANSYS at t=0.01 (c)
σ y computed by the IPG-DDA at t=0.02 (d)
σ y computed by ANSYS at t=0.02Figure 19
The zoomed deformations near the peak of the sliding plane solved out by the IPG-DDA. Four cases using different friction coefficients and time limits are illustrated. The inclination angleθ of the sliding plane is such that tan⁡θ=0.125.
(a)
tan ⁡ φ = 0 . 12 at t=0.01 s: slip expected (b)
tan ⁡ φ = 0 . 12 at t=0.02 s: slip expected (c)
tan ⁡ φ = 0 . 126 at t=0.01 s: stick expected (d)
tan ⁡ φ = 0 . 126 at t=0.02 s: stick expectedFigure 20
The contours of the block horizontal displacement componentsux (m) solved out by the IPG-DDA in the cases of using the different friction coefficients tan⁡φ=0.12<tan⁡θ and tan⁡φ=0.126>tan⁡θ at the time instants t=0.01 s and t=0.02 s.
(a)
tan ⁡ φ = 0 . 12 at t=0.01 s: slip expected (b)
tan ⁡ φ = 0 . 12 at t=0.02 s: slip expected (c)
tan ⁡ φ = 0 . 126 at t=0.01 s: stick expected (d)
tan ⁡ φ = 0 . 126 at t=0.02 s: stick expectedFigure 21
The contours of the block vertical displacement componentsuy (m) solved out by the IPG-DDA in the cases of using the different friction coefficients tan⁡φ=0.12<tan⁡θ and tan⁡φ=0.126>tan⁡θ at the time instants t=0.01 s and t=0.02 s.
(a)
tan ⁡ φ = 0 . 12 at t=0.01 s: slip expected (b)
tan ⁡ φ = 0 . 12 at t=0.02 s: slip expected (c)
tan ⁡ φ = 0 . 126 at t=0.01 s: stick expected (d)
tan ⁡ φ = 0 . 126 at t=0.02 s: stick expected


## 6. Conclusion

In this paper, we describe a new coupling scheme for modeling the contact effects of deformable blocks. The IPG method is introduced to improve the accuracy of the DDA through discretizing blocks into discontinuous elements, while the DDA is used to deal with element contacts. Two auxiliary criteria for determining element contact types are also developed. The coupling of the IPG method and the DDA takes place in the phase to assemble the global equilibrium equations. During this process, the IPG method contributes the submatrices derived from the variational formulations defined on elements, and the DDA provides the submatrices obtained by minimizing contact spring strain energies. Such mixed strategy leads to a unified system of equilibrium equations. In this way, both block deformations and contact displacements can be simultaneously solved out. Numerical experiments demonstrate that the proposed method is accurate and effective for contact simulations.From the computing formulations derived in this paper, we can find that our coupling scheme can return to the original DDA if no block is subdivided, and it can also become a standard DG-FEM when contacts are not considered. Therefore, the numerical examples falling in these two extreme cases are not conducted here, although we have used some of them to verify our computer program, including more friction examples which have been verified in plenty of literatures [43, 44].It is worth noting that there exist some restrictions of DDA contact model, such as the difficulties to choose appropriate contact spring stiffness and to find the convergent solution of open-close iterations. However, the proposed coupling scheme adopts the contact model inherited from the DDA rather than more accurate and stable linear complementarity formulation as reported in literatures [45–50] and the references therein. The reason lies in that the focus of this paper is the coupling scheme based on the original DDA but not its variants. So we will incorporate linear complementarity formulation into the DDA in our next work.



---

*Source: 6217679-2016-11-30.xml*