Mathematical Problems in Engineering

Mathematical Problems in Engineering / 2016 / Article

Research Article | Open Access

Volume 2016 |Article ID 6217679 |

Yue Sun, Xiangchu Feng, Jun Xiao, Ying Wang, "Discontinuous Deformation Analysis Coupling with Discontinuous Galerkin Finite Element Methods for Contact Simulations", Mathematical Problems in Engineering, vol. 2016, Article ID 6217679, 25 pages, 2016.

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

Academic Editor: Gerardo Severino
Received01 Aug 2016
Revised18 Oct 2016
Accepted19 Oct 2016
Published30 Nov 2016


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.

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 [58]. 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 [918]. 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 [2226], and coupling the DDA and finite element methods (FEMs) [2730].

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) [3638]. 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 [2226] 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 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.

2.1. Problem Definition

Assume the block with the boundary has the mass density . is subjected to the body force per unit area. consists of two mutually disjoint parts and , where the prescribed displacement and the surface traction are defined, respectively. Use to denote the displacement and for the Cauchy stress tensor of the block ; then the dynamics of the block at the time can be stated aswhere is the gradient operator and represents the unit outward normal vector of the boundary . Use to indicate the strain tensor, and there exist the following relationships:where 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. Let be a nondegenerate (the nondegeneracy means that there exists such that if 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 by , where and . As shown in Figure 1, we call as an interior edge if or an exterior edge if . Denote the edges of the element by . Then the finite element subspace on can be defined as [34]where denotes the space of square-integrable functions and is the space of the polynomials of the degree less than or equal to the positive integer . is full of discontinuous piecewise polynomials and is used as the test function space below.

2.3. DG-FEM Weak Form

According to the IPG method [34, 36, 37], the variational formulation of problem (1)–(3) is to find , such thatwhere the bilinear form :and the linear form :where and denote the penalty parameter and the shear modulus, respectively, is the length of the element edge , and denotes the unit outward normal vector to ; the average operator and the jump operator across are defined aswhere is the common edge of the elements and .

Substitute (8) and (9) into (7), and we havewhere , andwhere the operator “” indicates the dot product of two vectors.

3. Global Computing Scheme

Consider the -block system , . 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 , ,. Numerically discretize the displacement of the element bywhere denotes the unknown vector variable to be determined on the element . We adopt the following and :where denotes the barycenter coordinate of the element , , are the translational displacements at , is the rotation at , and , , denote the strain components of the element .

In every time step, the following global equilibrium equations are formed in a mixed strategy:where the stiffness matrix is partitioned by the submatrices , the right hand side consists of , and contains the unknowns to be determined on the th element, . 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 [2730]. 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 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].

3.1. Inertia Force

In the time step interval , substituting (19) into (13) leads towhere 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) aswhere and indicate Newmark parameters and and 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 and in every time step:where and refer to the finial velocity and acceleration, respectively. In the numerical examples of this article, we choose and 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 , and therefore the vector in (24) can be divided out from the equations. So we assemble on the right side of (24) to in the stiffness matrix of the global equilibrium equations (22), while and are added to 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 element by and , respectively. Substituting (21) into (4) and (5), respectively, leads towhere and for plane stress problems and and 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 havewhere is the area of the element . The matrix is assembled to the submatrix in the stiffness matrix of the global equilibrium equations (22).

3.3. Body Force

Assume the element bearing the body force . Substituting (19) into (16) giveswhere the integral on the right of the above equality is assembled to the right hand side of the global simultaneous equilibrium equations (22).

3.4. Interior Interface Stiffness

Assume the interior edge as the common boundary of the elements and . Rewrite in (15) into the summation below:whereUsing the notations of average and jump operators defined in (10) and (11), respectively, we can expand the expressions above into the formulations below:Then substitute (19) into (31), (32), and (33), respectively, and we havewhere the matrix is spanned by the unit outward normal vector to the edge :Actually, transfers the surface traction 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.

3.5. Displacement Boundary Constraint

Assume the exterior edge of the element is prescribed with the displacement boundary constraint . We can expand (15) as follows by using the notations of jump and average operators defined in (10) and (11):Substituting (19) into (38) giveswhere 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 towhere 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 on its exterior edge . Then, substituting (19) into (17) leads towhere the integral on the right is assembled to the submatrix 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 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 contacts with the edge . A normal contact spring with the stiffness and a shear contact spring with the stiffness are applied at the contact point . 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, . From Figure 3(b), we can find the fact that is the distance between and , and therefore we havewhere denotes the area of the triangle of , , and and indicates the length of the edge :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 during this time step: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 onto the edge , and therefore it can be determined as follows:Due to the fact and , , (52) becomesWhen 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 . Then substituting (21) into (53) giveswhere .

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:where , and denote the th component of and the th component of , respectively, and the matricesThe 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 Figure 4(a), assume that the vertex is sliding along the edge , and the normal contact spring with the stiffness is applied at the contact point . As shown in Figure 4(b), and are the sliding displacements of and along , respectively, and denotes the normal contact displacement. Then the normal contact force and the friction force arewhere is the sign function to determine the relative sliding direction, and is the friction angle. The friction force acts on both sides of the sliding surface. Consequently, the potential energy contributed by the friction force consists of the energy for and the energy for :Based on the fact that is the projection of the displacement of on the contact surface , while is the projection of the movement of on , we can compute and as follows:where is identical to the expression in (48) and can be approximated by defined in (50) when using small enough time step length. Then, substitute (19) and (21) into (60) and (61), respectively, and we haveNow, using the expressions of and in the above two equalities, we minimize with respect to and , respectively:where 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).

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 loop and the open-close iteration :(1)Initialize the time step counter ; then input data and set the maximum time step number .(2), 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].(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, and go to (); if false, modify the contact springs which violate the criterion, and reassemble the global equilibrium equations; then go to ().

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 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 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 Figure 1 as an example, where two blocks are in the surface-to-surface contact between the broken lines and . The letter is omitted in the figure. Here, and are located on the same node but belong to different elements, . Based on DDA contact detection criteria and (C1), elements contact between the edge couples: and , and , and , and and . These contacts consist of the vertex-to-edge contacts as enumerated in Table 1. Note that the vertex couples and , and , and , and , and , and and are not vertex-to-vertex contacts. Take and as an example to explain the reasons. According to (C2), is not allowed contact with the interior edge , while the situation that contacts with has been taken into account as the contact of the vertex with the edge .

Pair numberInvasion vertexEntrance lineInvasion vertexEntrance line


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 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 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 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 . 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 110.
Comparative Task 2. The results acquired from the IPG-DDA and ANSYS are compared on the horizontal displacement component of the right half of Block 1. Figures 7(a) and 7(b) draw the contour plots of computed by the two methods, respectively. It is shown that the displacement varying in the 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 of the right half of Block 1. The contour plots of 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 of Block 1.

Block 1Block 2

Plane stressYesYes
Mass density (Kg/m2)00
Young’s modulus (GPa)101000
Poisson’s ratio0.30


Element for Block 1PLANE82Constant stress element
Element for Block 2Pilot nodeDDA block
Contact elementCONTA172 and TARGE169
DOFs per element12/PLANE82 and 3/pilot6
Number of elements888 PLANE82s and 1 pilot888 elements and 1 block
Total DOFs106595334
Normal contact stiffness (N/m)114114
Shear contact stiffness (N/m)