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

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.

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][37][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][23][24][25][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.

Linear Elasticity Dynamics on a Single Block
Consider the block system containing 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.

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 Γ and Γ , 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 ∈ [0, ∞) can be stated asü 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: where C is a fourth-order elastic tensor, which is symmetric, positive definite, and point-wisely constant in Ω.

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. Let T ℎ = { 1 , 2 , . . . , Λ } be a nondegenerate (the nondegeneracy means that there exists > 0 such that if ℎ = diam( ) then contains a ball of radius ℎ in its interior. The operator diam( ) is the supremum of the distance between any two points in ) triangular subdivision of the block Ω. Denote the edges of all elements in T ℎ by E = { 1 , 2 , . . . , , +1 , . . . , }, where { 1 , 2 , . . . , } ⊂ Ω and Γ = { +1 , . . . , }. As shown in Figure 1, we call as an interior edge if 1 ⩽ ⩽ or an exterior edge if + 1 ⩽ ⩽ . Denote the edges of the element by . Then the finite element subspace on T ℎ can be defined as [34] D (T ℎ ) = {k ∈ 2 (Ω) : k| ∈ P ( ) , ∀ ∈ T ℎ } , (6) where 2 (Ω) denotes the space of square-integrable functions and P ( ) is the space of the polynomials of the degree less than or equal to the positive integer . D (T ℎ ) is full of discontinuous piecewise polynomials and is used as the test function space below.

DG-FEM Weak Form.
According to the IPG method [34,36,37], the variational formulation of problem (1) where the bilinear form (u, k) : and the linear form (k) : D (T ℎ ) → R: where and denote the penalty parameter and the shear modulus, respectively, | | is the length of the element edge ∈ E, and n denotes the unit outward normal vector to ; the average operator {⋅} and the jump operator ⟦ ⋅ ⟧ across are defined as where is the common edge of the elements + and − . Substitute (8) and (9) into (7), and we have where S ≡ (E \ Γ) ∪ Γ , and where the operator "⋅" indicates the dot product of two vectors.

Global Computing Scheme
Consider the -block system {Ω [1] , Ω [2] , . . . , Subdivide the blocks whose stresses are to be refined, and the rest intact blocks are regarded as block-type elements. Assume elements are created from the considering block system, including 1 , 2 ,. . . , . Numerically discretize the displacement u of the element by where { } denotes the unknown vector variable to be determined on the element . We adopt the following ( , ) and : where (̃,̃) denotes the barycenter coordinate of the element , 0 , V 0 are the translational displacements at (̃,̃),  In every time step, the following global equilibrium equations are formed in a mixed strategy: where the stiffness matrix is partitioned by the submatrices ∈ R 6 × 6 , the right hand side consists of ∈ R 6 × 1 , and ∈ R 6 × 1 contains the unknowns to be determined on the th element, , = 1, 2, . . . , . 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.
It 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][28][29][30]. The latter ones need to solve DDA and FEM equations alternatively, and in timedependent 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 entries and 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].

Inertia Force.
In the time step interval [ , + Δ ], substituting (19) into (13) leads to 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 { } in (23) as where and indicate Newmark parameters and { 0 } and { 0 } 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 Mathematical Problems in Engineering 5 recursive formulations to update { 0 } and { 0 } in every time step: where { 1 } and { 1 } 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 functions k , and therefore the vector { } T in (24) where = and ] = ] for plane stress problems and = /(1−] 2 ) and ] = ]/(1−]) for plane strain problems, in which the constants and ] indicate Young's modulus and Poisson's ratio of the element , respectively. The matrix [ ] above is known as the material matrix of . 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 where | | is the area of the element .

Body Force.
Assume the element bearing the body force b = ( , ) T . Substituting (19) into (16) gives where the integral on the right of the above equality is assembled to the right hand side { } of the global simultaneous equilibrium equations (22).

Interior Interface Stiffness.
Assume the interior edge as the common boundary of the elements and . Rewrite in (15) into the summation below: where Using the notations of average and jump operators defined in (10) and (11), respectively, we can expand the expressions above into the formulations below: 6

Mathematical Problems in Engineering
Then substitute (19) into (31), (32), and (33), respectively, and we have where the matrix [ ] is spanned by the unit outward normal vector ( , ) to the edge : Actually, [ ] transfers the surface traction ( , ) T to the combination of stress components in the following way: where denotes the stress tensor referring to (3). For the edge with the ending vertex coordinates ( , ) and ( , ), its outward unit normal vector can be computed by = − and = − .
Assemble the coefficient matrices in (34) into the stiffness matrix of the global equilibrium equations (22) according to the rules below: 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 and the integral transform technique.

Displacement Boundary Constraint.
Assume the exterior edge of the element is prescribed with the displacement boundary constraint u = ( , ) T . We can expand (15) as follows by using the notations of jump and average operators defined in (10) and (11): Mathematical Problems in Engineering 7 Substituting (19) into (38) gives where [ ] is the unit outward normal matrix introduced in (35). Assemble the coefficient matrices in (39) into the submatrix [ ] in the stiffness matrix of (22): Besides, substituting (19) into (18) leads to where the coefficient vectors are assembled to the right hand side { } of the global equilibrium equations (22); that is, 3.6. Surface Traction. Assume the element subjected to the traction t = ( , ) T on its exterior edge . Then, substituting (19) into (17) leads to where the integral on the right is assembled to the submatrix { } in the right hand side of (22).

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 vertexto-edge cases, but the first two can always be decomposed into vertex-to-edge contacts. As defined in Section 2.2, exterior element edges are portions of block boundaries, while interior elements edges are the others. We use the sign to indicate the set of the boundary of the element . As shown in Figure 3(a), assume that the vertex 1 ∈ contacts with the edge 2 3 ⊂ . A normal contact spring with the stiffness ⊥ and a shear contact spring with the stiffness ‖ are applied at the contact point 0 ∈ . Contact springs contribute the strain energies Π ⊥ and Π ‖ : where ℎ ⊥ and ℎ ‖ 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, ℎ ⊥ and ℎ ‖ have to be figured out in the first step. Assume the coordinates of the point at the beginning and the end of the current time step are ( , ) and ( , ), respectively, = 0, 1, 2, 3. From Figure 3(b), we can find the fact that ℎ ⊥ is the distance between 1 and 2 3 , and therefore we have where denotes the area of the triangle of 1 , 2 , and 3 and indicates the length of the edge 2 3 : If we choose small enough time step length, the last item on the right side of (47) can be ignored, and can also be approximated by the initial length of 2 3 during this time step: Mathematical Problems in Engineering Substitute (21), (49), and (50) into (46), and we have the following approximation of ℎ ⊥ : The shear contact displacement ℎ ‖ is actually the projection of the line 0 1 onto the edge 2 3 , and therefore it can be determined as follows: Due to the fact = − and V = − , = 0, 1, . . . , 3, (52) becomes When we use small enough time step length, the last item on the right side of the above equality can be ignored, and can be approximated by 0 . Then substituting (21) into (53) gives After substituting (51) into (44) and (54) into (45), we minimize the summation of Π ⊥ and Π ‖ with respect to { } and { } and have the following contact submatrices: Mathematical Problems in Engineering where , = 1, 2, . . . , 6, and denote the th component of and the th component of , respectively, and the 1 × 6 matrices 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).

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 Figure 4(a), assume that the vertex 1 ∈ is sliding along the edge 2 3 ⊂ , and the normal contact spring with the stiffness ⊥ is applied at the contact point 0 ∈ 2 3 . As shown in Figure 4(b), 0 and 1 are the sliding displacements of 0 and 1 along → 2 3 , respectively, and ℎ ⊥ denotes the normal contact displacement. Then the normal contact force and the friction forcêare where sgn(⋅) is the sign function to determine the relative sliding direction, and > 0 is the friction angle. The friction forcêacts on both sides of the sliding surface. Consequently, the potential energy Π contributed by the friction forceĉ onsists of the energŷ0 for and the energŷ1 for : Based on the fact that 0 is the projection of the displacement of 0 on the contact surface → 2 3 , while 1 is the projection of the movement of 1 on → 2 3 , we can compute 0 and 1 as follows: where is identical to the expression in (48) and can be approximated by 0 defined in (50) when using small enough time step length. Then, substitute (19) and (21) into (60) and (61), respectively, and we have Now, using the expressions of 0 and 1 in the above two equalities, we minimize Π with respect to { } and { }, respectively: where , = 1, 2, . . . , 6 and and denote the th component of and the th component of , respectively. The two vectors above are assembled to the right hand side of the global equilibrium equations (22). (1) Initialize the time step counter ← 0; then input data and set the maximum time step number * .

Implementation Aspects
(2) ← + 1, stop the program if > * , 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].  stresses, and velocities, and go to (2); if false, modify the contact springs which violate the criterion, and reassemble the global equilibrium equations; then go to (4).

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 vertexto-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 Figure 1. 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 openclose 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 Figure 1 as an example, where two blocks are in the surface-to-surface contact between the broken lines 1 , and + 4 5 and + 9 10 . These contacts consist of the vertex-to-edge contacts as enumerated in Table 1. Note that the vertex couples + 2 and − 7 , − 2 and + 7 , + 3 and − 8 , − 3 and + 8 , + 4 and − 9 , and − 4 and + 9 are not vertex-to-vertex contacts. Take + 2 and − 7 as an example to explain the reasons. According to (C2), + 2 is not allowed contact with the interior edge − 7 11 , while the situation that 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.

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 Figure 5(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 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 Figure 5(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        Table 4. Block 1 is subjected to the pressure 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 of the contact region is where is the radium of the cylinder and is Young's modulus; the normal contact pressure ( ) varying along the interface has the following analytical solution:  and the maximum contact pressure 0 takes place at the origin of the coordinate system: For Example 2, the theoretical half width of the contact region is ≈ 3.321276 mm, and the analytical maximum contact pressure 0 ≈ 383.103328 MPa. 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 Figures 9(b) and 9(c). The computational settings of the two methods are listed in Table 5. Besides, we adopt = 100 for the penalties used in (37a)-(37d), (40), and (42)    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 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 0 , 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 obtained by the IPG-DDA and ANSYS, while Figure 11 Figure 13(a) plots the contour distribution of the 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 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 solved out by the IPG-DDA and ANSYS together to compare them. Figure 14  of 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 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 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.
Example 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 Figure 15 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. Figure 15(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 = 1000 in (37a)-(37d), (40), and (42)        incline angle of the top surface of Block 2, and 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).  Figure 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. 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 -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. Figure 19 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)    IPG-DDA agree with the predictions based on kinetics very well.
In Figure 20, the distributions of the block horizontal displacement components solved out by the IPG-DDA are compared in the cases of using different friction coefficients.

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][46][47][48][49][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.