Advances in Materials Science and Engineering
Volume 2018 (2018), Article ID 6783791, 11 pages
https://doi.org/10.1155/2018/6783791

Research Article

A Simple Model for Elastic-Plastic Contact of Granular Geomaterials

Jian Wang,1,2,3 Qimin Li,⁴ Changwei Yang,¹ Yidan Huang,¹ and Caizhi Zhou⁵

¹School of Civil Engineering, Southwest Jiaotong University, Chengdu 610031, China
²MOE Key Laboratory of High-speed Railway Engineering, Southwest Jiaotong University, Chengdu 610031, China
³National Engineering Laboratory for Technology of Geological Disaster Prevention in Land Transportation, Southwest Jiaotong University, Chengdu, Sichuan 610031, China
⁴School of Water Resources and Environment, China University of Geoscience (Beijing), Beijing 100083, China
⁵Department of Mechanical Engineering, University of South Carolina, Columbia, SC 29208, USA

Correspondence should be addressed to Qimin Li; qiminli@163.com

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 as

where 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 as

where 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 as

where 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 × 10⁻⁴, and 1 × 10⁻³ m, respectively. And in the normal direction, element size of each zone changes gradually from 2 × 10⁻³ m to 2 × 10⁻⁴ m towards the contact area. The circumference is divided into 60 equal segments.

(a)

(b)

Figure 1: 3D finite element mesh for identical spheres. (a) 3D FE mesh for two identical spheres in contact. (b) Zone partition of the sphere.

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.

Table 1: Parameters used in the FEA analysis.

Number
(GPa)

(°)
(MPa)

1
10
0.10
20
5

2
20
0.20
30
10

3
30
0.30
40
20

4
—
—
—
30

5
—
—
—
40

6
—
—
—
50

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.

Figure 2: Contact pressure versus normal force.

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.

Figure 3: Initiation of plastic flow in spherical contact (a quarter of cross section from the center of the upper sphere). PEMAG, plastic strain magnitude.

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)

Figure 4: Yield surface for the linear D-P criterion versus V-M criterion and the linear D-P criterion versus M-C criterion: (a) D-P versus Mises criterion in the plane; (b) D-P versus M-C criterion in the aπ plane.

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.

Figure 5: Distribution of normal pressure on the contact surface for .

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.

Figure 6: Normal force-displacement curves for dimer stainless steel pairs.

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 as

where to are the regressed results for those coefficients individually, as shown in Table 2.

Table 2: Regression results for loading and unloading coefficients.

Parameters
a
B
m
n

0.5398
1.3087
2.3104
1.2377

1.0279E − 03
−8.7130E − 04
9.6205E − 03
−1.2352E − 03

−4.7213E − 02
4.4259E − 02
−2.8768E − 01
3.6574E − 02

7.6500E − 05
−7.6019E − 05
9.1244E − 04
−9.5000E − 05

−6.1476E − 03
6.4120E − 03
−7.1235E − 02
8.0222E − 03

5.5848E − 05
−3.8951E − 05
6.6345E − 04
−7.9491E − 05

−4.4023E − 03
3.3763E − 03
−4.8732E − 02
6.2875E − 03

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.

Table 3: Parameters used for interpolation.

Type
(GPa)

(°)
(MPa)
(N)
(mm)

1
10
0.15
25
15
9.06
0.00164

2
10
0.25
35
25
84.45
0.00706

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)

Figure 7: Interpolation of the loading coefficient a for E =10 GPa. (a) Interpolated sample space; (b) the predicted value.

(a)

(b)

Figure 8: Interpolation of the loading coefficient b for E = 10 GPa. (a) Interpolated sample space; (b) the predicted value.

(a)

(b)

Figure 9: Interpolation of the unloading coefficient m for E = 10 GPa. (a) Interpolated sample space; (b) the predicted value.

(a)

(b)

Figure 10: Interpolation of the unloading coefficient n for E = 10 GPa. (a) Interpolated sample space; (b) the predicted value.

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.

Figure 11: Normal force versus normal displacement for elastic-plastic deformation (type 1 in Table 3).

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.

Figure 12: Normal force versus normal displacement for elastic deformation (type 2 in Table 3).

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.

Table 4: Errors caused by comparing with FEA results.

Type
(%)
(%)
(%)
(%)
Method

1
1.92
0.50
9.31
0.75
Interpolation

1
4.05
0.66
34.30
1.47
Equation (19)

2
0.09
0.10
1.20
0.16
Interpolation

2
3.84
0.49
56.37
0.90
Equation (19)

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)

References

H. Kruggel-Emden, E. Simsek, S. Rickelt, S. Wirtz, and V. Scherer, “Review and extension of normal force models for the discrete element method,” Powder Technology, vol. 171, no. 3, pp. 157–173, 2007.

L. Vu-Quoc, X. Zhang, and L. Lesburg, “A normal force-displacement model for contacting spheres accounting for plastic deformation: force-driven formulation,” Journal of Applied Mechanics (ASME), vol. 67, no. 2, pp. 363–371, 2000.

M. R. Brake, “An analytical elastic-perfectly plastic contact model,” International Journal of Solids and Structures, vol. 49, no. 22, pp. 3129–3141, 2012.

D. Rathbone, M. Marigo, D. Dini et al., “An accurate force–displacement law for the modelling of elastic–plastic contacts in discrete element simulations,” Powder Technology,, vol. 282, pp. 2–9, 2015.

C. Thornton, S. J. Cummins, and P. W. Cleary, “An investigation of the comparative behaviour of alternative contact force models during inelastic collisions,” Powder Technology, vol. 233, pp. 30–46, 2013.

Z. Song and K. Komvopoulos, “An elastic–plastic analysis of spherical indentation: constitutive equations for single indentation unloading and development of plasticity due to repeated indentation,” Mechanics of Materials, vol. 76, pp. 93–101, 2014.

X. Zhang and L. Vu-Quoc, “An accurate elasto-plastic frictional tangential force-displacement model for granular-flow simulations: displacement-driven formulation,” Journal of Computational Physics, vol. 225, no. 1, pp. 730–752, 2007.

X. Zhang and L. Vu-Quoc, “A method to extract the mechanical properties of particles in collision based on a new elasto-plastic normal force-displacement model,” Mechanics of Materials, vol. 34, no. 12, pp. 779–794, 2002.

Abaqus/Standard, Version 2017. Dassault Systemes Simulia Corp, Abaqus/Standard, Providence, RI, USA, 2017.

W. F. Chen and G. Y. Baladi, Soil Plasticity: Theory and Implementation, Elsevier, New York, NY, USA, 1985.

L. Vu-Quoc and X. Zhang, “An elastoplastic contact force-displacement model in the normal direction: displacement-driven version,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 455, no. 1991, pp. 4013–4044, 1999.

L. Vu-Quoc, X. Zhang, and L. Lesburg, “Normal and tangential force-displacement relations for frictional elasto-plastic contact of spheres,” International Journal of Solids and Structures, vol. 38, no. 36-37, pp. 6455–6489, 2001.

Q. J Zheng, H. P Zhu, and A. B. Yu, “Finite element analysis of the contact forces between a viscoelastic sphere and rigid plane,” Power Technology, vol. 226, pp. 130–142, 2012.

E. Fjaer, R. M. Holt, P. Horsrud et al., Petroleum Related Rock Mechanics, Elsevier, Amsterdam, Netherlands, 2nd edition, 2008.

R. E. Goodman, Introduction to Rock Mechanics, Wiley, New York, NY, USA, 2nd edition, 1989.

S. Antonyuk, S. Heinrich, J. Tomas et al., “Energy absorption during compression and impact of dry elastic-plastic spherical granules,” Granular Matter, vol. 12, no. 1, pp. 15–47, 2010.

L. Kogut and I. Etsion, “Elastic–plastic contact analysis of a sphere and a rigid flat,” Journal of Applied Mechanics (ASME), vol. 69, no. 5, pp. 657–662, 2002.

L. R. Alejano and A. Bobet, “Drucker-Prager criterion,” Rock Mechanics and Rock Engineering, vol. 45, no. 6, pp. 995–999, 2012.

C. D. Wilson, “A critical reexamination of classical metal plasticity,” Journal of Applied Mechanics (ASME), vol. 69, no. 1, pp. 63–67, 2002.

E. Wang, T. On, and J. Lambros, “An experimental study of the dynamic elasto-plastic contact behavior of dimer metallic granules,” Experimental Mechanics, vol. 53, no. 5, pp. 883–892, 2013.

L. Y. Li, C. Y. Wu, and C. Thornton, “A theoretical model for the contact of elastoplastic bodies,” Proceedings of the Institution of Mechanical Engineers, Part C: Journal of Mechanical Engineering Science, vol. 216, no. 2, pp. 421–431, 2002.

I. Etsion, Y. Kligerman, and Y. Kadin, “Unloading of an elastic-plastic loaded spherical contact,” International Journal of Solids and Structures, vol. 42, no. 13, pp. 3716–3729, 2005.

R. Hemsley, in Interpolation on a Magnetic Field, Bristol University, Bristol, UK, 2009, Technical Report.