Advances in Materials Science and Engineering

Volume 2018, Article ID 6783791, 11 pages

https://doi.org/10.1155/2018/6783791

## A Simple Model for Elastic-Plastic Contact of Granular Geomaterials

^{1}School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China^{2}MOE Key Laboratory of High-speed Railway Engineering, Southwest Jiaotong University, Chengdu 610031, China^{3}National Engineering Laboratory for Technology of Geological Disaster Prevention in Land Transportation, Southwest Jiaotong University, Chengdu, Sichuan 610031, China^{4}School of Water Resources and Environment, China University of Geoscience (Beijing), Beijing 100083, China^{5}Department of Mechanical Engineering, University of South Carolina, Columbia, SC 29208, USA

Correspondence should be addressed to Qimin Li; moc.361@ilnimiq

Received 6 February 2018; Revised 13 April 2018; Accepted 3 May 2018; Published 3 June 2018

Academic Editor: Kaveh Edalati

Copyright © 2018 Jian Wang et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

We propose a simple elastic-plastic contact model by considering the interaction of two spheres in the normal direction, for use in discrete element method (DEM) simulations of geomaterials. This model has been developed by using the finite element method (FEM) and nonlinear fitting methods, in the form of power-law relation of the dimensionless normal force and displacement. Only four parameters are needed for each loading-unloading contact process between two spheres, which are relevant to material properties evaluated by FEM simulations. Within the given range of material properties, those four parameters can be quickly accessed by interpolating the data appended or by regression functions supplied. Instead of the Von Mises (V-M) yield criterion, the Drucker–Prager (D-P) criterion is used to describe the yield behavior of contacting spheres in this model. The D-P criterion takes the effects of confining pressure, the intermediate principal stress, and strain rate into consideration; thus, this model can be used for DEM simulation of geomaterials as well as other granular materials with pressure sensitivity.

#### 1. Introduction

DEM simulations play an indispensable role in exploring the multiscale relationship of particulate material properties, and it is also a feasible alternative to experimental investigations [1]. Interaction between two spheres, usually presented as force-displacement relationship, is the most fundamental problem in DEM simulations which dominates the motion of the particle system [2]. For dry, noncohesive granular media, viscoelastic models or hysteretic models are always used to describe their interparticle behavior, irrespective of the energy dissipation forms. The viscoelastic models are velocity damping dependent, while the hysteretic models are plasticity related. Because of the unphysical behavior of viscoelastic models, for example, the non-zero initial force at the beginning and end of the collision, hysteretic models have attracted more attentions in this field [3–6].

The existing hysteretic models were derived from either analytical methods or numerical simulations and designed for applications with regard to material properties as well as load level. Those models based on numerical results can describe the contact force-displacement relationship accurately during the elastic-plastic deformation; thus, they are termed “accurate models” [4, 7]. Also, they are empirical models as input parameters are fitting to material properties by regression. Many of those accurate models have some disadvantages as follows: (a) the complex implicitly equations are not easy to be implemented into DEM programs [1]; (b) it is difficult to get their input parameters [8]. Furthermore, the material parameters used in most accurate models are derived from the V-M yield criterion, which are not suitable for materials effected by confining pressure, intermediate principal stress, and strain rate, for example, geomaterials and other granular materials with pressure sensitivity.

Thus, this work proposes a simple and “accurate” model to describe the elastic-plastic behavior of granular geomaterials characterized by the D-P yield criterion. Also, a fast and accurate parameters-accessing method is provided based on the interpolation/regression method to make the model feasible to DEM simulations.

#### 2. Yield Criteria for Geomaterials

The pressure-dependent D-P yield criterion is used for geomaterials in this study and can be written aswhere is a pseudo–effective stress, is the slope angle of the linear yield surface in the stress plane (meridional plane), is the hydrostatic pressure, and is the effective cohesion of the material.

The flow potential, , for the linear Drucker–Prager model is defined aswhere is the dilation angle in the plane. Set , resulting in associated plastic flow, and assume that the ratio of the yield stress in triaxial tension to the yield stress in triaxial compression equal to 1, that is, the flow stress ratio [9], which implies that the yield surface of this model in the deviatoric principal stress plane is the V-M circle. Then, (1) can be rewritten aswhere is the first invariant of the stress tensor and is the Mises equivalent stress and defined by the second invariant of the deviatoric tensor as

By comparing with its two-parameter ( and ) form [10], we can obtain the following equation:

It can be deduced that

It should be noticed that the input friction angle in the linear Drucker–Prager model can be related to the friction angle from the triaxial test as follows:

If we nondimensionalize the vertical position in the sphere by the corresponding contact radius, that is, , then the principal stresses can be expressed as

Then, (4) can be rewritten as

By defining the amplitude function of the stress field of the D-P model as a three-parameter-dependent equation, we can get the following equation:

Then, we can get the initial yield point where is maximized with respect to by solving the partial derivative:

When we find the yield height, we can get the corresponding maximum contact pressure , and then the initial yield force and yield displacement can be expressed as (12) and (13), respectively:

#### 3. Finite Element Analyses

##### 3.1. Finite Element Model for Elastic-Plastic Contact

A 3D model in ABAQUS/Implicit is used to simulate the interaction between two identical spheres. As the contact is highly localized, with the ratio of contact radius to the sphere radius , only the region close to the contact area is considered in our calculations according to the Saint Venant principle as shown in Figure 1, which is similar to the model used in the work done by Vu-Quoc and coworkers [11, 12] and Zheng et al. [13]. The modeling domain was divided into three zones and adaptive mesh was used in our 3D model. The zones I and II were defined within the circular regions with the radius of 0.02*R* and 0.05*R*, respectively, and zone III contained the area outside the Zone II. The finite element mesh consists of 73100 8-node linear bricks (C3D8) and 76830 nodes. In radial direction, the element sizes of those three zones are 5 × 10^{−5}, 5 × 10^{−4}, and 1 × 10^{−3} m, respectively. And in the normal direction, element size of each zone changes gradually from 2 × 10^{−3} m to 2 × 10^{−4} m towards the contact area. The circumference is divided into 60 equal segments.