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.
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 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.
(a) Stacked three blocks and the implemented meshes
(b) Structure of the stiffness matrix
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–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 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.
(a)
(b)
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).
(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 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 .
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.
(a) Geometry sizes and the coordinate system
(b) Mesh for the right half of Block 1
(a) IPG-DDA result
(b) ANSYS result
(a) IPG-DDA result
(b) ANSYS result
Example 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 Figure 9(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 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 iswhere 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 takes place at the origin of the coordinate system:For Example 2, the theoretical half width of the contact region is , and the analytical maximum contact pressure .
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 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 with the coordinates . 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 , 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(b) draws the vertical displacement component 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 direction from the two methods, respectively. Figures 12(c) and 12(d) present the contour images of the displacement in the 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 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 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(a) draws the contour 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.
(a) Geometry sizes
(b) Triangular mesh
(c) Mesh detail nearby the contact region
(a) Horizontal displacement component
(b) Vertical displacement component
(a) computed by the IPG-DDA
(b) computed by ANSYS
(c) computed by the IPG-DDA
(d) computed by ANSYS
(a) computed by the IPG-DDA
(b) near the contact domain simulated by the IPG-DDA
(c) computed by ANSYS
(d) near the contact domain simulated by ANSYS
(a) computed by the IPG-DDA
(b) near the contact domain simulated by the IPG-DDA
(c) computed by ANSYS
(d) near the contact domain simulated by ANSYS
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(a), consider two linearly elastic trapezoidal blocks, Block 1 and Block 2. The computing time is limited in the range s. During the whole analysis, Block 1 is always subjected to a uniform pressure with the intensity N/m distributed on its top surface. At the initial time 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. 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 in (37a)–(37d), (40), and (42) for the IPG-DDA.
Comparative Task 1. Firstly, we compare the displacements of the point varying in the time interval obtained by the IPG-DDA and ANSYS. The point located at coordinates is one of the corners of Block 1. Figures 16(a) and 16(b) plot the computed displacement of in the and directions, respectively. In the two figures, the analytical solutions are computed based on Newton’s laws and the theory of kinematics as follows: where N/m is the traction intensity, is the top surface length of Block 1, /m2 is the mass density, m2 is the area of Block 1, is the 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).
Comparative Task 2. Next, we compare the distributions of the stress component obtained by the IPG-DDA and ANSYS. Figures 17(a) and 17(b) show the contour plots of at s figured out by the two methods, respectively. Figures 17(c) and 17(d) illustrate the distributions of at 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 solved out by the IPG-DDA and ANSYS are compared. Figures 18(a) and 18(b) give the contour distributions of at s obtained by the two methods, respectively. Figures 18(c) and 18(d) draw the contours of at 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 to the negative direction of the -axis. According to kinetics theory, a stick response is expected when the friction coefficient used in (58) is larger than . We run a couple of computations by using the IPG-DDA as follows: the first adopts the joint property with the friction coefficient such that a slipping response is predicted along the interface; and the second chooses the joint property with so that a stick state is expected.
Figure 19 compares the computed geometry configurations of the considered blocks with the friction coefficients and 0.126 by the IPG-DDA. Figures 19(a) and 19(b) plot the case with the friction coefficient at the time instants s and s, respectively, while Figures 19(c) and 19(d) illustrate the case with the friction coefficient at s and s, respectively. In the first case, we have , so a slipping response is expected. In the second case, there is ; 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 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. Figures 20(a) and 20(b) show the contours of when at s and s, respectively. In this case, the slipping response is expected owing to . 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 when at s and s, respectively. In this case, the stick state is predicted because of . 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.
Figure 21 compares the distributions of the block vertical displacement components obtained by the IPG-DDA for different friction coefficients. Figures 21(a) and 21(b) plot the contours of for at s and 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 for at s and s, respectively, and the stick state is expected under this fictional condition. Distinctly discontinuous distribution of 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 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.
(a) Problem geometry
(b) Mesh configuration
(a) The displacement component (m) along the -axis
(b) The displacement component (m) along the -axis
(a) computed by the IPG-DDA at
(b) computed by ANSYS at
(c) computed by the IPG-DDA at
(d) computed by ANSYS at
(a) computed by the IPG-DDA at
(b) computed by ANSYS at
(c) computed by the IPG-DDA at
(d) computed by ANSYS at
(a) at s: slip expected
(b) at s: slip expected
(c) at s: stick expected
(d) at s: stick expected
(a) at s: slip expected
(b) at s: slip expected
(c) at s: stick expected
(d) at s: stick expected
(a) at s: slip expected
(b) at s: slip expected
(c) at s: stick expected
(d) at 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.
Competing Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgments
This research is supported by the National Natural Science Foundation of China (Grant nos. 61271294 and 61471338).