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.02R and 0.05R, 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.
(a)
(b)
3.2. Material Properties
As sandstone and granite materials are two typical nonmetallic materials, the mechanical properties of those two materials listed in Table 1 were used in our simulations according to the work by Fjaer et al. [14]. According to Table 1, there are a total of 162 different combinations of values existing. To make reliable statistical analysis, all of 162 combinations were used as the input for our calculations, which generated sufficient samples to be regressed or interpolated.
3.3. Loading Level for Granular Geomaterials
To decide the loading level of normal force, sandstone is used for simulation, with material properties: . For more details of this type of material, refer to Goodman [15]. The applied normal forces increase to 200 N gradually, with the corresponding maximum contact pressure and maximum Mises stress shown in Figure 2. Assuming the contact pressure between spheres equals to the hydrostatic pressure resulting from the gravity of granular media, then the contact pressure under normal force of 100 N is equivalent to the maximum pressure generated by the granular aggregate (with density ) with a height of 1650 m. Thus, the loading level can be set to be 100 N, which is large enough to deal with real problems as well as extreme conditions. Considering the effect of confined stress field, the micro-yield strength is higher than the macroscopic yield stress, and that makes our load level safer [16]. In addition, the increasing maximum V-M stress recorded reflects the effect of strain hardening after initial yielding, and the maximum V-M stress is always lower than the corresponding maximum contact pressure.
With the increasing of normal force, the plastic zone expands at the core region of the sphere after yielding. When the plastic zone expands to the free surface, the constraint from the elastic materials disappears and unconstrained plastic flow forms. Based on the FEA calculations, the incipient of plastic flow can be identified through the evolution of plastic strain in the sphere. For the sphere of sandstone, plastic flow starts at the normal force approaching to 535 N (Figure 3), which is larger than the loading level of 100 N. It means that the inner elastic core surrounded by plastic deformation will not disappear in the end. Normally, geological granular materials do not exhibit fully plastic behavior during the deformation [17]; thus, only elastic and elastic-plastic deformation needs to be considered in our contact model.
3.4. D-P Criterion for Geomaterial Simulation
There are two advantages to use the D-P criterion in our model: (i) the D-P criterion can provide hydrostatic pressure-dependent yielding condition, and (ii) the D-P criterion can make the calculation converge smoothly.
The shear failure mechanism described by the von Mises criterion is independent of the stress level (i.e., the hydrostatic pressure), which is mainly used to define the plastic deformation of metallic material with constant yield strength and can hardly apply to geomaterials, as shown in Figure 4(a).
(a)
(b)
For nonmetallic materials, the D-P criterion provides a smooth yield surface compared with the Mohr–Coulomb yield surface, the cross section of which in a deviatoric plane (π-plane) is an irregular hexagon with sharp corners that may cause convergence problems for numerical simulation. Although the Drucker–Prager criterion can provide a simple and smooth yield surface, it may not predict the failure accurately when one or more principal stresses are tensile [18]. It is worth to mention that the issue raised by tensile principal stresses in the Drucker–Prager criterion will not appear in our calculations, as the primary condition to use the D-P criterion is naturally satisfied due to the two dry, cohesionless spheres.
Based on the previous work done by Alejano and Bobet [18], both Mohr–Coulomb criterion and Drucker–Prager criterion can accurately predict the failure of nonmetallic materials for triaxial experiments around the triaxial compression stresses, that is, . If the value of differs from , then the Mohr–Coulomb criterion will underestimate the strength of the geomaterials and the Drucker–Prager criterion will overestimate it with the increase of intermediate principal stress. Actually, the intermediate and the minor principal stresses are always identical, that is, , for the spherical contact studied in this work, so those two criteria are automatically satisfied in our calculations.
To validate the Drucker–Prager criterion, the responses of those three models with equivalent yield strength under the same loading level have been analyzed and compared.
Following (6), the relationship between yield strength () used for the V-M model and uniaxial compressive strength used for the D-P model can be expressed as
The cohesion property used for the Mohr–Coulomb model is related to as follows:
The material parameters for those three criteria used for comparison are set to be as follows:(a)For the D-P criterion, (b)For the M-C criterion, , cohesion , and tension cutoff (c)For the V-M criterion, .
As far as the maximum contact pressure and contact area were concerned, the distributions of normal traction predicted by the Drucker–Prager criterion and the Mohr–Coulomb criterion are nearly the same, except the slight difference between their shapes in Figure 5. It is also found that the influence of flow stress ration K on the results of the Drucker–Prager criterion is negligible when comparing the result of and . However, the result predicted by the von Mises criterion deviated a lot from the other two curves which match well with the real pressure distribution within the contact area.
According to the above analysis, the Drucker–Prager criterion is the best one to analyze the interaction between nonmetallic spheres, while those results from the von Mises criterion can hardly transplant to nonmetallic materials.
3.5. Model Calibration
To validate to the D-P yield criterion used in the particles’ contact model, comparisons are made between the numerical and experimental results in Figure 6. It is worth to mention that the D-P criterion is not only applicable to most of the geomaterials but also suitable to some metallic materials exhibiting strain hardening behaviors, such as 2024-T351 aluminum alloy [19] and stainless steel 302 [20]. Except loading level and degree of plastic deformation, there is no fundamental difference between particles made of geomaterials and strain-hardened metals, as far as the contact force-displacement curve is concerned. The force-displacement curves of dimer bead pairs made of stainless steel 302 with diameters of 9.53 mm and 12.7 mm tested by the Split-Hopkinson pressure bar in low velocity are used for calibration. With the D-P yield criterion adopted, the values of input parameters for numerical model and material properties are taken from the experiment results [20]. Overall, the calculated results matched well with the experimental results. Moreover, the influence of dimensions of spheres on the test results can be eliminated largely by nondimensionalizing the normal displacement by the initial yield displacement , especially at the stage without obvious plastic flow. It means that the size effect of spheres in contact can be avoided by nondimensionalization for the same type of material. Accordingly, a nondimensionalized contact model can be developed based on numerical results to show the general relationship of normal force and displacement.
4. Results and Discussion
4.1. A New Normal Contact Model
According to the analysis in Section 3.3, the contact spheres of granular geomaterials hardly undergo fully plastic deformation, thus only elastic and elastic-plastic deformation needs to be considered. In addition, there is no pop-in behavior at the transition from the elastic regime to the elastic-plastic regime in the force-displacement curves. Thus, a single continuous function can be used to describe the force-displacement curve at the loading stage instead of a piecewise approach developed by other models [1, 21]. Relationships of nondimensionalized force and displacement are explored to eliminate the influence of size effect and to make a symmetric pattern for the equation. We perform a nonlinear least-squares fitting of the loading and unloading data for the elastic-plastic case. It is found that power-law functions are the best approximation to disrobe the force-displacement curves in the dimensionless form comparing with other relationships, such as exponential, polynomial, and parabola relationships, for their coefficients of determination (i.e., R-square) approaching 1 and the root mean squared errors (RMSEs) keep the smallest. All those 162 combinations listed in Table 1 have been tested, and the fitted parameters are appended in Supplementary Materials (available here).
Dimensionless force-displacement curves at the loading stage can be expressed with power-law relation as follows:where and are the contact load and normal displacement at the inception of plastic deformation, respectively, and and are the loading coefficients to be determined.
Inspired by Etsion et al. [22] and Song and Komvopoulos [6], force-displacement curves at the unloading stage can also be expressed as the similar form as loading curves:where is the residual displacement after unloading and and are the unloading coefficients.
The residual displacement can be found from the initial unloading process; that is, it can be predicted from the beginning of the unloading process. During the loading stage, the normal displacement of sphere center is gradually increased up to a desired maximum , and the contact force reaches its maximum . Then, the unloading process is initialized, and the residual displacement can be deduced from (17) as follows:
4.2. Regression Method for the Parameters
It is found that all of the four parameters have linear or nonlinear influence on the loading and unloading coefficients separately. The nonlinear effects of those parameters on the targeted coefficients can be approximated by quadratic polynomial fitting. To prepare for the regression of a multilinear model, the quadratic term was linearized and treated as a newly separate variable. After that, the stepwise regression has been performed by using the Matlab package, which can add or remove terms from a multilinear model based on their statistical significance. The empirical relationship between material properties and loading and unloading coefficients can be expressed aswhere to are the regressed results for those coefficients individually, as shown in Table 2.
R-squares for those four fitted functions are 0.7068, 0.7665, 0.5166, and 0.7517, respectively.
4.3. Interpolation Method for the Parameters
Both the Kriging and natural neighbor interpolations can be used to predict the loading and unloading parameters accurately. More details about those interpolation methods can be found in Hemsley [23]. 3D, 4D, or 5D interpolation may be employed according to the parameters selected in this data space, that is, the space formed by . The interpolation process can be described in the following three steps: Step 1: identify the minimum interval of interpolating points and discretize the data space accordingly Step 2: interpolate values to the discretized points by the proper method Step 3: search the predicted value at the target point locating on the grid of interpolation space.
To validate the new model (16)–(18), two types of materials with different mechanical properties listed in Table 3 have been chosen for interpolation by means of the natural neighbor interpolation method, which has been proved to be exact and it is continuous everywhere except at sample points [23]. Contacting spheres () with material as type 1 suffered plastic deformation, but it would probably stay in the elastic regime with material as type 2 under the maximum normal force of 100 N.
Figures 7–10 present the 4D interpolated results and corresponding slice for searching the target value (type 1 in Table 3) of those four needed parameters.
(a)
(b)
(a)
(b)
(a)
(b)
(a)
(b)
There is a wider span for the unloading coefficient in Figure 7 than the loading coefficient a in Figure 9, which means the unloading coefficient changes more acutely and irregularly than the latter, with the ratio of the variance to the minimum value 238% versus 29% for E = 10 GPa. Therefore, it is harder to capture the variance of m by the linearized regressing method, which may be responsible for the low R-square in (19).
Both the loading coefficient and the unloading coefficient are below 1.5 and in the similar range, as shown in Figures 8 and 10.
4.4. Comparing Fitted Results with FEM Results
With the model shown in Figure 1, materials in Table 3 and normal loading level of 100 N, normal force-displacement curves gained from FEM and predicted by the proposed model are put together in Figures 10 and 11. We can see that FEA results are reliable for the testified model.
As shown in Figure 12, the loading-unloading curve predicted by the interpolation method closely agree with the FEA results for elastic-plastic deformation, which means that the coefficients of the model proposed can be directly obtained from interpolation of the data samples supplied without doing time-consuming and complicated finite element analysis. However, the curve originating from the regression function can be treated as a coarse approximation to FEA results when there is lack of additional information. It also shows that the curve of elastic materials is always steeper than that of elastic-plastic materials as we expected.
The normal force-displacement curves for the Hertz theory, FEA, and interpolation are totally overlapped for elastic deformation, which is an excellent evidence of the validation of the interpolation method used in Figure 10. There is a little deviation of the regressive results to analytical and FEA results even for elastic material, which should be adjusted for actual use.
Coefficient errors caused by different data-processing methods by comparing with FEA results are shown in Table 4, with the materials used in Table 3. Two main trends can be found: (a) the errors caused by the regression function are larger than those caused by the interpolation method and (b) the errors of elastic-plastic materials are larger compared with elastic materials by use of the interpolation method. It also shows that the main source of errors for the regression method originates from the unloading coefficient , as analyzed in Section 4.3.
5. Conclusions
By taking the advantages of the Drucker–Prager yield criterion, a new elastic-plastic contact model has been developed for nonmetallic materials based on FEA results. This new model presented the following characters:(a)Considering the absence of plastic flow for granular geomaterials and no pop-in transition from the elastic regime to the elastic-plastic regime, power-law relationships of dimensionless force and displacement can be described by nonlinear least-squares fitting for the loading and unloading process between spheres in contact. They are concise and easy to be implemented to DEM codes.(b)Four parameters of the model are relevant to material properties taken from FEM results, which can be quickly accessed by interpolation or regressing equations instead of simulation or any other assumptions. The application and accuracy of the model can be expanded and enhanced with the enlargement of parameters database appended.(c)The D-P criterion is used to describe the yield behavior of contacting spheres in the model, which takes the effect of confining pressure, the intermediate principal stress, and the strain rate into consideration; thus, this model can be used for DEM simulation of geomaterials as well as other granular materials with pressure sensitivity.
Data Availability
The data used to support the findings of this study are included within the article and supplementary materials.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
The authors appreciate the financial support from the National Natural Science Foundation of China (no. 51308474) and the Fundamental Research Funds for the Central Universities (no. 2682017CX005).
Supplementary Materials
A data pool named “DataPoolForInterpolationMethod.xls” is provided with 162 samples calculated by FEA with the form of a table. It is used for interpolating the coefficients of loading and unloading models (i.e., (14) and (15)) in 4D space of parameters E, v, φ, and σc. For the material in the range of E = 10–30 GPa, = 0.10–0.30, φ = 20–40°, and σc = 5–50 MPa, each coefficient (i.e., a, b, m, and n) can be quickly accessed directly or by Kriging or natural neighbor interpolation. The yield strength as well as contact radius at the yield moment is also provided for each sample, which makes it more convenient for manipulation. (Supplementary Materials)