Repeated loading model for elastic–plastic contact of geomaterial

Abstract

A new nonlinear hysteretic model with considering the loading, unloading, and reloading processes is developed based on Drucker–Prager yield criterion and finite-element analysis. This model can be used for multiple repeated elastic–plastic normal direction contact problems between two identical spherical geomaterials. After examining the influence of material properties, strain hardening, and loading histories, we found that the hysteretic phenomena (represented by residual displacement and plastic work) become weak after the first cycle, and the subsequent cycles step into elastic shakedown state eventually. A critical number of cycles can be used to estimate the state of ratchetting, plastic shakedown, as well as elastic shakedown. It also found that the subsequent curves will be stiffer than the previous ones, especially when the yield strength is high and ratchetting effect is not strong. This new model can be used for a wide range of geomaterials under different loading levels, and it can also be extended to describe the constitutive behavior of spheres under earthquake as well as aftershocks.

Introduction

Interaction between two spheres is the most fundamental problems for discrete geological material simulation, among which repeated loading and unloading is especially important for dynamic action, such as compacted discrete rock mass or soil under seismic force.

Usually, very simple linear models of contact, like linear viscoelastic models, are used in Discrete Element Methods, which account for more than 60% of all cases. ¹ The main disadvantage of non-hysteretic models is their artificial assumption of force–displacement behavior, which would not guarantee their accuracy to describe the overall macroscope response of particle systems under repeated loading in elastic–plastic regime.

With regard to the interaction without cohesion, it has been realized that hysteretic models consisted of springs and sliders would be more physically intuitive than other types of models with dampers for the energy dissipation. The contact zone is under locally plastic deformation, rather than the artificial damper behavior. ²

Although many contact models of the hysteretic type have been developed for corresponding materials separately,^3–7 relatively few attentions have been paid to the unloading and reloading responses of elastic–plastic deformation especially, which is essential for cyclic loading simulation between spheres.

The main difficulties for cyclic loading model with hysteretic law result from the definition of stiffness coefficients and assessment of the evolution of stress state during the unloading and reloading processes. Literatures concerned include the elastic–plastic indentation of rigid sphere to compliant base and interference between two flexible spheres, which could be mutually emulated substantially.

Generally speaking, the understanding of unloading process goes through elastic assumption to elastic–plastic verification. Assuming the unloading behavior purely elastic, Johnson ⁸ proposed the analytical model of unloading for elastic–plastic indentation. Later, the model had been modified by Li and Gu ⁹ based on the superposition of Steuermann solutions instead of the original contact surface curvature. That modification made it more closely fit with the numerical results with large plastic deformation. Yan and Li ¹⁰ studied the deformation of a rigid sphere under cyclic indentation. In their study, an elastic–perfectly plastic half-space was adopted and they found that the deformation during unloading from a contact beyond the elastic limit was perfectly elastic and the reloading curve was the same as the unloading curve. The material of the half-space used by them represented a typical steel material, with high elastic modulus and yield stress. However, it is still unclear whether the unloading and reloading response would be different for other materials with low elastic modulus and yield stress.

Song and Komvopoulos ¹¹ examined the unloading behavior of indented elastic–plastic materials as well as the evolution of plasticity due to four consecutive indentation cycles by finite-element analysis (FEA) and found the occurrence of plastic deformation during unloading was attributed to the effect of the residual stresses in the elastic material surrounded by the plastic regime. Also, they deduced that elastic shakedown behavior (i.e. dynamic creep, creep associated with fatigue, or deflection stability) would be characterized by the elastic–plastic deformation for material with high effective elastic modulus-to-yield strength ratio (E*/Y) under repeated loading–unloading cycle, while it would keep ratcheting deformation for material with low E*/Y.

Concerning the spherical interference, Etsion et al. ¹² studied the unloading response of an elastic–plastic sphere contacting with a rigid plane by FEA and noticed that the purely elastic unloading was just one of the extreme case, and they concluded that all subsequent load–unload process after the initial loading–unloading cycle would be fully elastic. They proposed an unloading model of the power law relation with nonconstant exponent, which was also been adopted by Song and Komvopoulos. ¹¹ Jackson et al. ¹³ found high residual stresses remained after an elastic–plastic hemispherical contact being unloaded. Kadin et al. ¹⁴ presented the numerical results of two successive loading and unloading cycles on the same model as Etsion et al. ¹² and inferred that the secondary plastic flow might reach a steady state and the successive loading–unloading cycles become fully elastic from the evolution of equivalent von Mises stress distribution. Later, Zait et al. ¹⁵ extended the unloading contact model of Etsion et al. ¹² including stick contact condition.

All the above-mentioned results are mainly fit for metallic material following the von Mises yield criterion (i.e. Maxwell–Huber–Hencky–von Mises theory), which are not suitable for more complicated geo materials with hydrostatic pressure dependence. Besides, the reloading model based on reliable FEA results has not been put forward, which is crucial for cyclic loading problems.

Although Sadd et al. ¹⁶ developed a nonlinear hysteretic contact law including loading, unloading, and reloading stages for dry cohesionless granular media simulation based on photo-elastic experiment, the inter-particle contact behavior did not exceed the elastic regime substantially.

Kostek ¹⁷ presented a nonlinear mathematical model describing the hysteresis of dry contact of rough planar surfaces loaded in the normal direction, which is hard to transform to spherical contact.

Therefore, the main purpose of this article is to develop unloading and reloading models for geomaterial under cyclic loading and with elastic–plastic deformation. Also, the influence of material properties and loading histories on the evolution of elastic–plastic deformation will be examined. This article is structured as follows. FEA model and yield criterion are described in sections “Numerical model” and “Material constitutive model,” respectively. Schematic contact model is presented in section “Contact model,” and loading model of the first cycle, unloading model, and reloading model are discussed in sections “Loading model,”“Unloading model,” and “Reloading model” individually. Influence of material properties and loading history on loading–unloading curves is discussed in sections “Influence of material properties” and “Influence of loading history,” respectively. Finally, conclusions are drawn in section “Conclusion.”

Contact models

Contact model

Figure 1 presents the contact force F and displacement δ curve for elastic–plastic contact between two spheres. The loading process O → A can be predicted accurately by hysteretic model $F=k_n{\delta }^{1.5}$ , with the nonlinear spring stiffness k_n changing with different loading level, and k_n can be expressed as piecewise functions of displacement $\delta$ , that is, $k_n=f(\delta )$ at elastic, mixed elastic–plastic, and fully plastic regimes or just one regression function for the whole domain (equations (8)–(10)). A full description of the loading process of elastic–plastic contact is given in section “Loading model.”

Figure 1.

Schematic illustration of contact force versus displacement of loading, unloading, and reloading processes of elastic–plastic contact.

The contact loads in the unloading process A → B and reloading process B → C are affected by the initial loading level $F_{\max }$ , ${\delta }_{\max }$ , and residual displacement level ${\delta }_{\mathrm{res}}$ , so the unloading model and reloading model can be written in a similar form, such as equations (11) and (15).

The schematic model assumes that the reloading behavior is elastic. The reloading response is elastic shakedown under repeated loading, which can be simplified to elastic behavior when the incremental plastic strain and plastic work of the subsequent cycles are negligible, according to the shakedown theorems. This assumption will be further illuminated in section “Influence of loading history” in the following.

The schematic model presented is similar to the model proposed by Sadd et al. ¹⁶ in terms of the schematic form, while the latter is dependent on the previous loading and unloading paths. In the presented reloading model, only the final states of the previous loading and unloading processes, that is, the maximum loading level and the residual displacement of the previous unloading, are concerned, which makes it simpler without losing accuracy.

Numerical model

The interaction of two identical spheres (ball–ball interaction) is equivalent to a sphere acting on a rigid flat surface (ball–wall interaction) due to the symmetry of ball–ball interaction, while the displacement of ball–wall interaction is half of ball–ball interaction. In order to improve the computing efficiency without losing the accuracy, an axisymmetric ball–wall model is adopted from the nonlinear finite-element package ABAQUS 6.14.

The finite-element mesh shown in Figure 2(a) and (b) represents an axisymmetric model of hemisphere (with radius R = 0.1 m) with 32,828 four-node bilinear quadrilateral elements (CAX4R element). The hemisphere is divided into three different mesh density zones. Zone I, within 0.1R from the sphere tip, has as small mesh size as 0.0001 times the yield contact radius $a_c$ to achieve fine discretization. Zone II lays outside the 0.1R distance from the sphere tip vertically, but within 0.6R distance from the sphere center horizontally. The mesh size increases monotonically from the bottom to the top, as well as from the center to the edge of the sphere. Thus, Zone III has the coarser mesh than the other two parts. Besides, the rigid plane is modeled with a one-dimensional two-node axisymmetric rigid link (RAX2 element), with the size gradually increasing from the interaction point to the far end.

Figure 2.

(a) Schematic model of spherical contact and (b) finite element meshing.

The boundary conditions are also shown in Figure 2(a). The rigid plate is fixed on the bottom, thus only compliant deformation of the sphere allowed. The axisymmetrical boundary is applied to the center of the hemisphere, which prevents the moving of nodes in the radial direction. Contact force is applied to the center of the hemisphere. The vertical displacements of those nodes at the top are constrained by multi-point constraints (MPCs), and the interaction between ball and the wall is frictionless and set to be “surface-to-surface” discretization. The normal contact force F increases from 0 N to the target value with a constant rate in 3 min and then decreases back to 0 N steadily in 1 min, which makes a loading–unloading cycle.

The contact model has been verified by comparing the calculated results with the analytical output of Hertz solution in the elastic regime, as well as the analytical solution of the incipient yield force and displacement from the D-P criterion. Data for sandstone have been used for those verifications. Basically, FEA results closely agree with the Hertzian solution in the elastic regime, and the difference between the numerical and analytical results of the elastic–plastic constitutive model at the incipient yield point is less than 5%, thus the contact model presented is suitable for further exploration in the elastic–plastic regime, according to Vu-Quoc et al. ⁴

The material properties of Table 1 can be referred for the uniaxial compression test of sandstone, as shown by Goodman. ¹⁸

Table 1.

Model parameters.

Parameters	Values
Elastic modulus, E	10 GPa
Poisson’s ratio, υ	0.38
Internal frictional angle, φ	$27.8^{\circ }$
Uniaxial compression strength, σc	10 MPa

Material constitutive model

Yield criterion and hardening rule

The chosen yield criterion is the Drucker–Prager (D-P) linear yield function with associated plastic flow, which can be written as

$$f=q-p\tan \beta -d=\sqrt{3J_2}+\frac{I_1}{3}\tan \beta -d=0$$

(1)

where q is the generalized shear stress and defined by the second invariant of the deviatoric tensor J₂ as $q=\sqrt{3J_2}$ p is the mean stress, β is the slope angle of the linear yield surface in the p–q stress plane, and I₁ is the first invariant of the stress tensor. d defines the position of the initial yielding surface and the development of subsequent yielding surfaces, as shown in Figure 3. For the associated plastic flow, the slope of the linear yield surface equals to the dilation angle $\psi$ in the p–q plane, that is, $\psi =\beta$ .

Figure 3.

Schematic of hardening and flow for the linear D-P model in the p-q plane.

If hardening is defined by the uniaxial compression test, then the initial yield surface can be defined as

$$d=(1-\frac{1}{3}\tan \beta ){\sigma }_c$$

(2)

where ${\sigma }_c$ is the uniaxial compression yield stress.

The strain-hardening curve as well as its parameters for geomaterial can be obtained and implemented to predict stress–strain responses closely agreeing with the experimental tests, according to the method proposed by Jiang and Wu. ¹⁹ The actual exponential strain-hardening relationship can be idealized and approximated by bilinear model. ²⁰ Thus, the linear isotropic hardening curve can be shown as Figure 4, with the ratio of hardening modulus to elastic modulus m. m is also called strain-hardening parameter. ²¹ For the perfectly plastic case, m equals to 0, and for the elastic case, m equals to 1. For simplification and universality, only the fully plastic case has been chosen for the contact model proposed, but the influence of strain-hardening parameter on contact response will be discussed.

Figure 4.

Stress σ versus strain ε for elastic linear-strain-hardening plastic model.

As shown in Figure 4, the yield stress with isotropic hardening effect for the elastic linear-strain-hardening plastic material can be written as

$$\sigma {\prime }_c={\sigma }_c+\frac{\mathrm{mE}}{1-m}{\epsilon }_p$$

(3)

where E is the elastic modulus, ${\sigma }_p$ is the plastic stress, and ${\epsilon }_p$ is the plastic strain.

Thus, the subsequent yielding surface in Figure 3 can be defined as

$$\mathrm{d}\prime =(1-\frac{1}{3}\tan \beta )({\sigma }_c+\frac{\mathrm{mE}}{1-m}{\epsilon }_p)$$

(4)

Flow rule

The associated flow potential, g, for the linear D-P model can be written as

$$g=q-p\tan \psi$$

(5)

Isotropic strain-hardening defined in uniaxial compression can then be modeled in terms of the incremental plastic strain as

$$\mathrm{d}{\epsilon }_{\mathrm{ij}}^p=d{\overline{\epsilon }}_p\frac{1}{(1-\frac{1}{3}\tan \psi )}\frac{\partial g}{\partial {\sigma }_{\mathrm{ij}}}$$

(6)

where ${\overline{\epsilon }}_p$ is the equivalent plastic strain; a scalar measure of plastic strain tensor defined through time integration of the plastic strain rate ${\stackrel{·}{\epsilon }}_{\mathrm{ij}}^p$ is represented as

$${\overline{\epsilon }}_p=\underset{0}{\overset{t}{\int }}{\left(\frac{2}{3}{\stackrel{·}{\epsilon }}_{\mathrm{ij}}^p{\stackrel{·}{\epsilon }}_{\mathrm{ij}}^p\right)}^{1/2}\mathrm{dt}$$

(7)

Results and discussion

The loading, unloading, and reloading processes of spherical contact have been studied, considering the influence of material properties and loading history.

Spherical contact model

Loading model

The contact model for geomaterial can be summarized as equations (8) and (9)

$$F\left(\frac{\delta }{{\delta }_y}\right)=\{\begin{array}{c}\hfill k_e{\left(\frac{\delta }{{\delta }_y}\right)}^{\frac{3}{2}},0\le \frac{\delta }{{\delta }_y}<1\hfill \\ \multicolumn{1}{c}{k_{\mathrm{ep}}{\left(\frac{\delta }{{\delta }_y}\right)}^{\frac{3}{2}}+{\Delta }_{\mathrm{ep}},1\le \frac{\delta }{{\delta }_y}<\frac{{\delta }_p}{{\delta }_y}}\\ \multicolumn{1}{c}{k_p{\left(\frac{\delta }{{\delta }_y}\right)}^{\frac{3}{2}}+{\Delta }_p,\frac{\delta }{{\delta }_y}\ge \frac{{\delta }_p}{{\delta }_y}}\end{array}$$

(8)

where stiffness coefficient is

$$\{\begin{array}{c}\hfill k_e=\frac{4}{3}{E}^{*}{R}^{{*}^{\frac{1}{2}}}{\delta }_y^{\frac{3}{2}}\hfill \\ \multicolumn{1}{c}{k_{\mathrm{ep}}=\frac{p_1{\left(\frac{\delta }{{\delta }_y}\right)}^2+p_2\left(\frac{\delta }{{\delta }_y}\right)+p_3}{{\left(\frac{\delta }{{\delta }_y}\right)}^2+q_1\left(\frac{\delta }{{\delta }_y}\right)+q_2}}\\ \multicolumn{1}{c}{k_p=m\left(\frac{\delta }{{\delta }_y}\right)+n}\end{array}$$

(9)

In this example, the pseudo stiffness coefficients in these three regimes can be regressed and expressed as equation (10), with the compensation of errors ${\Delta }_{\mathrm{ep}}=0.2185$ and ${\Delta }_p=4.6549$ , respectively

$$\{\begin{array}{c}\hfill k_e=6.9274\hfill \\ \multicolumn{1}{c}{k_{\mathrm{ep}}=\frac{5.273\times {\left(\frac{\delta }{{\delta }_y}\right)}^2+158.3\times \left(\frac{\delta }{{\delta }_y}\right)+42.61}{{\left(\frac{\delta }{{\delta }_y}\right)}^2+26.74\times \left(\frac{\delta }{{\delta }_y}\right)+8.304}}\\ \multicolumn{1}{c}{k_p=-0.004326\times \left(\frac{\delta }{{\delta }_y}\right)+5.676}\end{array}$$

(10)

The compliance curve predicted by this model is shown in Figure 5. The present model adopts the Hertzian solution in the elastic regime, as well as the best-fit curves in the elastic–plastic regime and in the fully plastic regime, with a small tuning at the boundary of their transition to eliminate the deviation caused by finite sliding contact-tracking approaches. The present model closely agrees with the FEA results. Moreover, it is more physically intuitive than those regressive models totally dependent on FEA results or analytical models constrained by some assumptions for the plastic deformation.

Figure 5.

The curve for the proposed model compared to the FEA results.

Unloading model

The numerical results as well as the empirical relation of unloading process initiated from elastic–plastic deformation (contact force F lies between the incipient yield force F_y and the incipient fully yield force F_p) and fully plastic deformation (contact force F lies above the incipient fully yield force F_p) are shown in Figure 6. As claimed by Song and Komvopoulos, ¹¹ plastic deformation also occurs during unloading, which makes the load–displacement behavior of unloading process nonlinear and strongly affected by the initial loading level $F_{\max }$ (i.e. the maximum loading force) and displacement level ${\delta }_{max}$ (i.e. the maximum displacement). Thus, the normalized contact force F/F_max during unloading is dependent on ${\delta }_{max}$ and residual displacement ${\delta }_{\mathrm{res}}$ . The relationship between unloading normal force F and normal displacement $\delta$ can be written as non-dimensional form as follows

$$\frac{F}{F_{max}}={\left(\frac{\delta -{\delta }_{\mathrm{res}}}{{\delta }_{max}-{\delta }_{\mathrm{res}}}\right)}^n$$

(11)

The exponent n can be expressed as the function of ${\delta }_{\max }$ in elastic–plastic regime and fully plastic regime (see equation (12))

$$n=-0.04846{\left(\frac{{\delta }_{\max }}{{\delta }_y}\right)}^{0.3091}+1.569$$

(12)

where ${\delta }_y$ is the incipient yield displacement for the material dealt with here, and it can be obtained following the method proposed by Vu-Quoc et al. ⁴ when replacing the von Mises law with D-P criterion. Actually, the incipient yield displacement for the material with D-P criterion not only depends on the Poisson’s ratio as that of von Mises criterion but also the angle of inner friction.

Figure 6.

Contact force versus displacement of unloading process (F is the contact force, F_y is the incipient yield force, and F_p is the incipient fully yield force).

It should be noticed that n = 1.521 when ${\delta }_{\max }={\delta }_y$ , which is a little bit higher than the Hertz solution (n = 1.5) for purely elastic loading because of the computational error caused by the contact-tracking approaches used for large deformation. Similar phenomena have also been found by Zhao et al., ²² Kogut and Etsion, ²³ and Yan and Li. ¹⁰ The best fitted coefficients used here can be compensated or adjusted slightly to make it agree with the analytical solution better. The error analysis and compensation between FEA and analytical results, as well as the analysis of incipient yield displacement for geomaterial, is a separate topic which is outside the scope of this article.

The residual displacement can also be predicted by the level of previous unloading displacement according to FEA results (see Figure 7(a) and equation (13))

$$\frac{{\delta }_{\mathrm{res}}}{{\delta }_{max}}=0.08536({\left(\frac{{\delta }_{max}}{{\delta }_y}\right)}^{0.3597}-1)$$

(13)

Figure 7.

(a) Normalized residual displacement and (b) normalized plastic work versus normalized maximum displacement of unloading process.

Similarly, the power function of normalized maximum displacement used to express the normalized plastic work closely agrees with the finite-element method (FEM) results (see Figure 7(b))

$$\frac{W_p}{W_{\mathrm{total}}}=0.07844({\left(\frac{{\delta }_{max}}{{\delta }_y}\right)}^{0.3644}-1)$$

(14)

where W_p is the plastic work and W_total is the total work the normal force has done.

Reloading model

Figure 8(a) presents the normalized contact force–normalized displacement results for reloading process with the maximum initial loading (F_max) is 6000 N, which attempts to be higher than the subsequent loading level. With the increase in the maximum initial loading, the residual contact area increases, which makes the reloading curve hardened.

Figure 8.

(a) Normalized contact force versus normalized displacement and (b) contact force versus displacement of reloading process (solid line: FEM results; Dashed line: fitted by equation (12)).

Similar to the unloading model, the reloading model can also be expressed as a power function as the unloading model

$$\frac{F}{F_{max}}={\left(\frac{\delta -{\delta }_{\mathrm{res}}}{{\delta }_{max}-{\delta }_{\mathrm{res}}}\right)}^m$$

(15)

The exponent m can be expressed as the function of ${\delta }_{max}$ in elastic–plastic regime and fully plastic regime, which decreases with the increase in previous loading level (see equation (16))

$$m=-0.106{\left(\frac{{\delta }_{max}}{{\delta }_y}\right)}^{0.2928}+1.688$$

(16)

The fitted function closely agrees with the FEM results, especially when the contact load is not too large, as shown in Figure 8(b).

It should be noticed that the reloading process predicted by equations (8), (11), and (13) is only suitable for those interaction without seriously plastic flow or particle breakage. Parameters in equations (10), (12), and (14) are related to material properties, which can be determined from experimental calibration tests or numerical calculations.

Influence of material properties

We explored the influence of elastic modulus E, compression yield stress ${\sigma }_c$ , inner frictional angle $\varphi$ , and strain-hardening parameter m on the unloading behavior. Figure 9 compares the normalized loading–unloading curve from low to high elastic modulus. We can see that the higher the elastic modulus, the wider area enveloped by the loading–unloading curve.

Figure 9.

Normalized contact force (the ratio of normal force F to maximum force F_max) versus normalized displacement (the ratio of displacement δ to maximum normal displacement δ_max) of loading and unloading processes with elastic modulus in the range of 1–100 GPa for F_max = 6000 N.

Figure 10 compares the results from different normalized parameters. We can see that the nonlinear spring stiffness monotonically increases with the increase in elastic modulus and decreases with the increase in inner frictional angle and strain-hardening parameter, but almost keep the same with the changing of compression yield stress. Normally, the incipient yield force F_y only depends on the inner frictional angle and the compression yield stress. Thus, it is reasonable to see that the maximum normalized load remains constant for different elastic modulus and strain-hardening parameter, while vary with the inner frictional angle and compression yield stress in Figure 10. Similar to the findings from Shankar and Arjun, ²⁴ our numerical results of isotropic hardening are consistent with analytical illustration in equation (3) and Figure 4, that is, the material behavior changes from fully plastic to perfect elastic with the strain-hardening parameter m transforming from 0 to 1.

Figure 10.

Normalized contact force (the ratio of normal force F to incipient yield force F_y) versus normalized displacement (the ratio of normal displacement δ to incipient normal displacement δ_y) of loading and unloading processes (a) with elastic modulus from 1 to 100 GPa, (b) with inner frictional angle from $5^{\circ }$ to ${30}^{\circ }$ , (c) compression yield stress from 5 to 20 MPa, and (d) strain-hardening parameter from 0 to 0.8 for F_max = 6000 N.

It should be noticed that the phenomenon of overlapping of the force–displacement curves has not been seen here as the metallic material following von Mises criterion with the change in material properties, which just indicates that the model derived from D-P criterion is more sensible to material properties than that from von Mises criterion. Actually, the whole overlapping of the force–displacement should never been expected for two different kinds of material, even for metallic material.^2,11 According to many unloading and reloading simulating results, the unloading and reloading models are still validated for other geomaterial with different frictional angle and yield stress, after the constants in equations (12) and (16) being adjusted by numerical or test results, which is the specific property of geomaterial with D-P criterion and can be easily implemented.

Influence of loading history

To get full insight into the effect of loading history on elastic–plastic response of geological material, 10 consecutive cycles are simulated for maximum loading force F_max = 3000 N (repeated loading). Seven other consecutive cycles with the maximum loading force monotonically increasing from 100 to 3000 N (monoincreasing loading) and decreasing from 3000 to 100 N separately have also been calculated (monodecreasing loading). Considering the sensitivity of material properties on residual displacement and plastic work, two levels of uniaxial compression stress are chosen for analysis (i.e. σ_c = 10 and 1 MPa), instead of the effective elastic modulus used by Song and Komvopoulos. ¹¹

Repeated loading

The relationship between contact force and displacement of the first six cycles are shown in Figure 11(a) and (b). When the uniaxial compression yield stress is at a high level (e.g. 10 MPa), the subsequent reloading and unloading behaviors are hardened after the initial loading process, which is similar to the hardening effects shown in Figure 10(d), but the residual displacement continuously increases with monodecreasing increment, as shown in Figure 11(a).

Figure 11.

Contact force versus displacement for repeated loading: (a) σ_c = 10 MPa, (b) σ_c = 1 MPa, and (c) time history of loading.

According to the shakedown theorems,^20,25,26 when both of the incremental plastic strain $\delta {\epsilon }_p$ and incremental plastic work $\delta W_p$ are nonzero, ratcheting occurs, so whether the value of the incremental plastic strain or incremental plastic work is zero can be used to identify the state of plastic response. When the uniaxial compression yield stress is at a lower level (e.g. 1 MPa), the subsequent repeated cycles lead to ratcheting behavior, which makes the incremental plastic strains repetitively accumulated (see Figure 11(b)).

Because of the rapid decrease in residual displacement and plastic work after the initial cycle, the ratios of residual displacement and plastic work of subsequent cycles to that of the first cycle are examined for clearness (see Figure 12(a) and (b)). When the compression yield stress of the sphere is 10 MPa, those two types of ratio are power law relation of the cycle number N separately, which are given by equations (17) and (18). When the compression yield stress is 1 MPa, those two types of ratio are fitted by equations (16) and (17) individually

$$\frac{{\delta }_{\mathrm{res}}}{{\delta }_{\mathrm{res}}^1}=23.73N^{-1.83}-0.3112$$

(17)

$$\frac{W_p}{W_p^1}=37.57N^{-1.663}-0.7669$$

(18)

where ${\delta }_{\mathrm{res}}^1$ and $W_p^1$ are residual displacement and plastic work of the initial cycle, respectively.

Figure 12.

(a) Ratio of subsequent residual displacement and (b) plastic work to those of the first cycle separately versus cycle number for repeated loading.

The critical cycles, that is, the cycle when the incremental plastic strain or incremental plastic work is zero, can be deduced by those ratio functions. Both critical cycles calculated by equations (17) and (18) are 11, while the critical cycles calculated by equations (19) and (20) are 17 and 23, respectively. Following the shakedown theorems, the sphere with higher compression yield force more easily transforms to elastic shakedown behavior. For the sphere with lower compression yield force, it will experience ratchetting (i.e. $\delta {\epsilon }_p\ne 0$ , $\delta W_p\ne 0$ ) in the first few repeated cycles (N < 17) and then turns into plastic shakedown state (i.e. $\delta {\epsilon }_p=0$ , $\delta W_p\ne 0$ , 17 ≤ N < 23) and it will end with elastic shakedown state finally (i.e. $\delta {\epsilon }_p=0$ , $\delta W_p=0$ ) after tens of cycles (N ≥ 23). The repeated action with tens of cycles between spheres with high stress level is not common, so ratcheting plays a dominant role in the plastic deformation for most of the case

$$\frac{{\delta }_{\mathrm{res}}}{{\delta }_{\mathrm{res}}^1}=3.165N^{-0.4468}-0.9059$$

(19)

$$\frac{W_p}{W_p^1}=6.59N^{-0.341}-2.288$$

(20)

It should be noticed that equations (17)–(20) are fitted from results of FEA, which show the influence of material properties on shakedown behavior of spherical contact. No limitation has been seen for the power law relationship, but the critical cycle will change with different yield strength.

The response of spherical contact starts from plastic flow and ends with elastic shakedown, in agreement with the finding of Kadin et al. ¹⁴ The detailed paths from plastic deformation to elastic shakedown response have been revealed, which may reflect the dominant effect of shakedown or ratchetting, as is shown in Figure 11(a) and (b).

Also, the geomaterial with D-P criterion used here does not so easily transform to elastic shakedown state under repeated loading as the metallic materials with von Mises criterion used by Kadin et al. ¹⁴ and Song and Komvopoulos. ¹¹ For their results, the earliest elastic shakedown behavior emerges at the second cycle, whereas in this study more than 11 cycles are needed.

Monotonically decreasing loading

The loads applied for monotonically decreasing loading are as follows: after the initial cycle of maximum loading force 3000 N, the maximum reloading forces from cycle 2 to cycle 7 are 2500, 2000, 1500, 1000, 500, and 100 N.

The reloading responses for materials with high or low level of compression yield strength follow the same way: the reloading process is rapidly stabilized to elastic shakedown response, as shown in Figure 13. For most materials exhibiting strain-hardening behavior, it is usual to see the stiffer reloading curve. The strain hardening defines that the strength of materials increases with the strain. As we adopted the strain-hardening equation in our model, the numerical effect is not applied to get the stiffer reloading curve in our calculations, which is caused by accumulated plastic deformation generated from previous loadings in the materials.

Figure 13.

Contact force versus displacement for monodecreasing loading: (a) σ_c = 10 MPa, (b) σ_c = 1 MPa, and (c) time history of loading.

Figure 14(a) and (b) further validates such elastic shakedown behavior. The incremental plastic work as well as incremental residual displacement reduces to zero quickly in 2–5 reloading cycles, which means the emergence of elastic shakedown response. It also shows that the material with lower compression yield strength more easily transforms to elastic shakedown state. Moreover, it should be noticed that the plastic work and residual deformation caused by reloading process are negligible to those of the initial loading process (the former are 3.1% and 2.4% of the latter for σ_c = 10 MPa, as well as 0.26% and 0.09% for σ_c = 1 MPa, separately; see Figure 14). Thus, the reloading process can be simplified to elastic response reasonably, and that is the base of our proposed reloading model. This finding validates the postulation of nonlinear hysteretic normal contact law presented in Figure 1, in agreement with that of Sadd et al. ²⁷

Figure 14.

(a) Ratio of plastic work and (b) residual displacement to those of the first cycle separately versus cycle number for monodecreasing loading.

The results of monodecreasing loading can be extended to a more general situation because of the elastic shakedown properties. If the initial loading level is higher than the subsequent loading level, no matter monodecreasing or not, the reloading process can be approximately simplified as elastic behavior, that is, the extended monodecreasing loading.

Monotonically increasing loading

The monotonically increasing loading process with maximum load from 100 to 3000 N is shown in Figure 15. Due to the effect of residual plastic deformation of the previous loading, the consecutive loading curve is hardened, and the stiffness of the virtual nonlinear spring is strengthened, comparing with the single load of the same level of maximum load. Also, the residual displacement of each cycle is a little bit larger than single load, while the plastic work caused by each cycle is obviously lower than single load in our example (see Figure 16(a) and (b)). It should be noticed that the plastic deformation continuously expands with the subsequent loading, and no shakedown behavior is expected for mono-increased loading process.

Figure 15.

(a) Contact force versus displacement for mono-increased loading for σ_c = 10 MPa and (b) time history of loading.

Figure 16.

(a) Normalized residual displacement and (b) normalized plastic work versus normalized displacement for σ_c = 10 MPa.

Seismic loading

The seismic action is cyclic or repeated loading essentially and can be decomposed as the superposition of single loading–unloading process without considering the difference between transient process and steady state, so it is meaningful to examine the relationship among the natural dynamic actions and those three types of idealized loading histories above. With the absolute values of maximum amplitude being scaled to the same value (i.e. 3 m/s ² ), the acceleration time histories of three famous earthquake records are shown in Figure 17(a)–(c). The three earthquake events referred are El Centro earthquake with 6.9 M_w occurred on 18 May 1940, Kobe earthquake (Great Hanshin earthquake) with 6.9 M_w happened on 17 January 1995, and Wenchuan earthquake with 7.9 M_w took place on 12 May 2008. With regard to the wave shape, the maximum acceleration amplitude of both El Centro earthquake and Kobe earthquake occurs quickly after the first arrival of P waves (less than 7% and 19% of the whole time recorded, respectively), so the effect of seismic action before the maximum amplitude is small and thus they can be represented as the extended monodecreasing loading process mentioned above. As for the Wenchuan earthquake, the maximum absolute value of acceleration happens a little later (about 25% of the whole time recorded), which shows the limited hardening effect as repeated loading or monoincreasing loading. Moreover, because of the elastic shakedown properties of monotonically decreasing loading mentioned above, the unloading and reloading models presented are also effective for the action of aftershocks, which are not negligible for their magnitude and amount.

Figure 17.

Seismograms of (a) El Centro earthquake, (b) Kobe earthquake, and (c) Wenchuan earthquake (the time of maximum amplitude is marked by dashed line).

Conclusion

The loading, unloading, and reloading behavior of spherical contact for elastic–plastic geomaterial has been analyzed by FEA calculations. The influence of material properties and loading history on the unloading and reloading responses has been examined. Based on numerical results and theoretical analysis, the main conclusions can be drawn as follows:

1. Hysteretic phenomena (represented by residual displacement and plastic work) decrease quickly after the first cycle, and the subsequent cycles will step into elastic shakedown state eventually, thus the hysteretic effect of subsequent cycles can be omitted reasonably.

2. A critical number of cycles can be used to estimate the state of ratchetting, plastic shakedown as well as elastic shakedown when repeated loadings are high.

3. The subsequent curves will be stiffer than the previous ones, especially when the yield strength is high and ratchetting effect is not strong.

4. This model presented can be extended to describe the deformation behavior of spheres under earthquake as well as aftershocks by considering the impact effect.