Numerical and experimental investigation of an elastoplastic contact model for spherical discrete elements

A contact model for the normal interaction between elastoplastic spherical discrete elements has been investigated in the present paper. The Walton–Braun model with linear loading and unloading has been revisited. The main objectives of the research have been to validate the applicability of the linear loading and unloading models and estimate the loading and unloading stiffness parameters. The investigation has combined experimental tests and finite element simulations. Both experimental and numerical results have proved that the interaction between the spheres subjected to a contact pressure inducing a plastic deformation can be approximated by a linear relationship in quite a large range of elastoplastic deformation. Similarly, the linear model has been shown to be suitable for the unloading. It has been demonstrated that the Storåkers model provides a good evaluation of the loading stiffness for the elastoplastic contact and the unloading stiffness can be assumed as varying linearly with the deformation of the contacting spheres. The unloading stiffness can be expressed in a convenient way as a function of the Young’s modulus and certain scaling factor dependent on the dimensionless parameter defining the level of the sphere deformation.


Introduction
Nowadays, the discrete element method (DEM) is used for modelling various particulate and nonparticulate materials such as soils, rocks, ceramics. The material in the DEM is represented by a large assembly of particles interacting among one another with contact forces. The particles can be of arbitrary shape, but spherical particles are often chosen because of the simplicity and computational efficiency of the numerical algorithm. The spherical particles will be considered in the present study.
The normal contact between particles in the DEM is very often modelled assuming an elastic-type force-displacement relationship. A simple contact model consisting of a linear spring in parallel with a viscous damping element proposed in the pioneering work by Cundall and Strack [4] is still very popular in DEM simulations [3,7,16]. A nonlinear elastic interaction is included in the contact models based on the Hertz theory [5,11,25].
In many applications, however, particle deformation due to contact cannot be treated as purely elastic. Because of the contact force concentration, yielding at the contact zone between two spheres made from ductile materials, e.g. metals, may occur at a relatively low loading [17]. In such cases, a partial irreversibility of interparticle penetration should be included in the contact model in the DEM. Various elastoplastic models have been proposed for the DEM, e.g. [10,15,[20][21][22][23][24].
The model proposed by Thornton [22] considers both the elastic and plastic ranges of the contact interactions. Initially, the contact is considered assuming the Hertzian elastic model. The loading is switched to the plastic regime when the yield criterion is satisfied. The yield criterion used by Thornton [22] is based on the assumption that the contact stress distribution after yielding can be obtained by "cutting off" the Hertzian stress distribution. The Thornton's model is used in different applications which require consideration of plastic effects in the interparticle interaction, e.g. [2]. It has been, however, found out that the Thornton's model underestimates significantly the contact force in the range of plastic loading, cf. [15,23]. Improved elastoplastic contact models analogous to that of Thornton have been proposed by Rathbone et al. [15] and by Vu-Quoc et al. [23]. In both cases, the improvements were based on the results of the finite element analyses of the contacting spheres.
The range of elastic loading is often very small in comparison with subsequent plastic loading, cf. [17]. Then, it can be neglected and the contact deformation can be treated as plastic from the contact initiation. This is consistent with the rigid-plastic or rigid-viscoplastic material model. A model assuming rigid-plastic behaviour according to the Hollomon stress-strain-hardening curve has been developed by Storåkers et al. [20,21]. The Storåkers contact model has been successfully used for discrete element modelling of powder compaction [12].
An issue of great importance in the elastoplastic contact models is the choice of a suitable unloading model. The Hertzian elastic unloading with a modified contact curvature proposed by Thornton [22] is usually assumed in the elastoplastic contact models, e.g. [15,23].
As it has been shown by Pasha et al. [14], a rigorous treatment of the elastoplastic contact deformation is often not necessary to represent behaviour of the granular material, therefore simple linear elastoplastic models such as that proposed by Walton and Braun [24] are still important for many practical applications. The elastoplastic loading and elastic unloading in the Walton-Braun model are governed by linear force-displacement relationships with different slopes, which ensures a residual overlap of the particles when the contact force drops to zero. No tensile forces are allowed in the Walton-Braun model. The elastoplastic model developed by Luding [10] can be considered as an extension of the Walton-Braun model to the adhesive contact.
Despite extensive research and apparent simplicity of the problem, the question of applicability ranges and accuracy of linear elastoplastic contact models in the DEM is not fully answered. The present work is aimed at numerical and experimental investigation of validity of the linear Walton-Braun-type elastoplastic model and suitability of analytical formulae for evaluation of model parameters defining the loading and unloading stiffness.
Discrete element modelling of powder metallurgy processes, including compaction, sintering and cooling is a practical motivation of the present work, therefore the formula proposed by Storåkers for powder compaction [8,13] will be verified for evaluation of the loading stiffness. A special interest will be paid to the unloading behaviour. The linear approximation of the unloading relationship will be checked against numerical and experimental data. The linear contact model in the DEM is consistent with linear macroscopic relationships, e.g. [9]. A linear unloading model is searched as a suitable microscopic counterpart of the linear thermoelastic macroscopic description of the unloading of the sintered material during cooling.
The outline of the present paper is as follows. In Sect. 2, a brief formulation of the contact problem in the discrete element method is given, and the formulation of the considered elastoplastic model: the Walton-Braun model with the linear loading and unloading is presented. The Storåkers model proposed for modelling the loading and the linear model for the elastic unloading are introduced. Section 3 presents laboratory tests performed to provide experimental data for validation of the analytical contact models and for determination of material properties. Compression of steel balls between two parallel flat plates has been carried out to obtain the force-displacement relationships for a large range of elastoplastic deformations. Uniaxial compression tests for special nonstandard samples have been performed to obtain stress-strain curve for finite element method simulations. Section 4 describes finite element method simulations of the problem analysed experimentally. Validation of the analytical contact model for the elastoplastic loading is shown in Sect. 5 using the experimental and numerical results. Modelling of the elastic unloading is considered in detail in Sect. 6.

Problem formulation
The normal contact interaction between two spherical particles i and j with radii R i and R j ( Fig. 1) is considered. The particles' positions are defined with the position vectors of Fig. 1 Definition of the two spheres contact their centroids, x i and x j . The contact force exerted on the i-th particle by the j-th is denoted by F i j , and by the Newton's third law the contact force exerted on the j-th particle F ji satisfies the relation: Further on, we take F = F i j = −F ji as the force representing the contact interaction between the particles. The contact force vector F can be projected on the unit normal vector n: where the unit normal vector n is assumed to be directed radially inwards to the centre of the particle i: We employ the so-called soft-contact DEM formulation, in which particles are treated as pseudo-rigid elements and the deformation is assumed to be localized at the contact zone. The effect of the local deformation is obtained by allowing a small overlap of the particles h ( Fig. 1): The amount of the overlap h is used for the determination of the contact force. Please note that the convention adopted for the vector n and the overlap h is such that the contact compressive force F is positive and so is the overlap h when the particles are in contact.

Walton-Braun model
The model proposed by Walton and Braun [24] assumes a linear force-overlap relationship, but the unloading slope (stiffness) is higher than the loading slope (stiffness), which leads to a certain residual irreversible overlap when the force drops to zero. This allows us to treat this model as an elastoplastic one. The force as a function of the particle overlap is plotted in Fig. 2. The contact force is given by: where h 0 is the residual irreversible overlap, and k L and k U are the loading and unloading stiffness, respectively, which should satisfy the condition: The residual overlap h 0 representing the plastic deformation of the contacting particles can be easily obtained as The reloding path follows the unloading path until the maximum overlap is achieved and the loading path is reactivated. The unloading stiffness k U can be taken as constant or variable described by a linear function of the maximum force, F max , or the maximum overlap, h max , achieved before unloading, so that: (8) or where S and B are certain constants.

Elastoplastic loading model
The elastoplastic loading in the presented above Walton-Braun model will be described by a special (linear) case of the Storåkers model [20,21], which has been proposed for a general viscoplasticity combining strain-hardening plasticity and creep. A simplified formulation without strain rate effects assumes rigid-plastic properties of the i-th and j-th particle materials which are defined by the Hollomon stress-strain relationships where σ (a) 0 , a = i, j, are material constants and m is the hardening exponent. The normal interaction force F is given by the following equation: where R * is the effective radius defined in terms of the particle radii, R i and R j h is the particle overlap given by Eq. (4), the parameter c 2 is given by [6]: and the equivalent Hollomon constant σ * 0 is defined as follows For the ideal plasticity (without hardening) and spheres made of the same material, when σ (a) = σ Y and m = 0, Eq. (11) is reduced to, cf. [13]: where c = 1.43. Comparing Eqs. (5) and (15), the loading stiffness k L can be written as:

Elastic unloading model
The Storåkers model has been derived neglecting elastic deformation, cf. [8]. In such a model, the unloading would be governed by the rigid behaviour (no change of deformation during the unloading). Assuming that the loading curve is valid for an elastoplastic material, the Storåkers model has been combined with the elastic unloading according to the Hertz model [12]. Here, we will assume the elastic unloading according to the linear relationship The stiffness k U will be evaluated assuming that the spring modelling contact elasticity at unloading is equivalent to an elastic bar of a nonuniform cross-sectional area (Fig. 3), consisting of two segments, with the lengths and the cross-sectional areas where 0 ≤ α i , α j ≤ 1 is the coefficient defining the area of the segment as a fraction of the particle cross-sectional area. The system of the two bar segments can be treated as two springs connected in series with the stiffnesses given by where E i and E j are the Young's moduli of the materials of the segments (or of the particles) i and j.
The equivalent stiffness k U of the system of two springs is given by, cf. [16]: Introducing the relationships (20) and (21) into the formula (22) and assuming E i = E j = E and α i = α j = α, we obtain the expression for the equivalent stiffness k U in the following form: where R * is the effective radius defined by Eq. (12). For equal size particles (R i = R j = R), Eq. (23): The purpose of further studies will be to determine the coefficient α suitable for evaluation of the unloading stiffness k U .
subjected to a contact pressure inducing a plastic deformation of the spheres. Two kinds of tests have been carried out. Firstly, compression of steel balls between two parallel flat plates has been performed in order to obtain experimental data on the contact forces during loading, unloading and reloading. Then, uniaxial compression tests have been carried out for the samples cut out from the balls in order to determine the stress-strain curve for finite element simulations.

Compression of steel balls
Forged balls made of steel grade DC04 with diameters of 30 mm, 40 mm, and 50 mm shown in Fig. 4 have been compressed between two flat plates. The experimental set-up is shown in Fig. 5. It can be easily noticed that this problem is equivalent to that of two equal balls subjected to contact under compressive loading (see Fig. 6) provided that the plates are ideally rigid. In both cases, the respective geometric quantities are related by the same equation: which allows us to derive force-displacement relationships valid for both cases. The balls after the tests with visible flattened contact areas are shown in Fig. 7. The force-displacement curves for the loading and unloading for the balls of different diameters are plotted in Fig. 8. It can be seen that in the investigated range of loading with a good accuracy, the curves can be approximated with straight line segments. It has been shown by Rojek et    al. [17] that the size-dependent force-displacement relationships for the contact between two equal balls can be easily transformed to the size-independent relationships between the pressure-type and strain-type parameters, F/4R 2 and h/2R, respectively. The force-displacement curves scaled in this way are plotted in Fig. 9. It can be observed that the scaled curves for the balls of different sizes merge pretty well. This shows that the scaled relationships can be used for different particle sizes provided that the material properties are the same. Some of the tests have been performed carrying out several loading-unloading-reloading (l-u-r) cycles. The scaled force-displacement curve for a ball of 30 mm diameter subjected to four l-u-r cycles is plotted in Fig. 10. It can be observed that the reloading curves coincide ideally with the

Determination of stress-strain curves
Nonstandard small cylindrical samples (of 5 mm diameter and 8 mm height) have been cut out from one of balls. Figure 11a shows the ball, the cylinder cut out from the ball, and one of the samples before and after compression test. The set-up of the compression test is shown in Fig. 11b. The stress-strain curve obtained from the compression tests is plotted in Fig. 12.

Finite element method simulation
The contact of an elastoplastic sphere against a rigid plane surface has been simulated using the finite element method. The analysed problem can be considered as quasistatic and it can be solved using either static or dynamic formulation. The explicit dynamic FE framework which is commonly used for metal forming processes [18,19] has been chosen here. The isotropy of material properties has been assumed. The geometry, loading and boundary conditions are symmetric about an axis of rotation; therefore, the problem can be solved using two-dimensional finite elements keeping the features of the three-dimensional description. A 2D axisymmetric model has been prepared for the ABAQUS/Explicit FE program [1]. A hemisphere with a radius of 15 mm has been discretized nonuniformly with 3175 four-node bilinear axisymmetric quadrilateral CAX4R finite elements with the reduced integration and hourglass control. The finite element model is shown in Fig. 13. A magnified detail in Fig. 13b illustrates a mesh refinement in the contact area where high stress and strain gradients are expected. Frictionless contact conditions have been assumed between the hemisphere and the plane. Compressive loading has been introduced by applying the vertical displacement to the top plane surface of the hemisphere without restrain- ing the horizontal motion. Similarly, the unloading has been controlled kinematically. In order to ensure quasistatic conditions, the loading and unloading have been applied using sinusoidal velocity profile with a sufficiently low-velocity amplitude. Maximum velocity 0.1 m/s has been taken for the loading and 0.02 m/s for the unloading.
Standard values of the Young's modulus and Poisson's for a steel have been taken, E = 2 · 10 5 MPa and ν = 0.3, respectively. The isotropic strain-hardening plasticity has been assumed with the stress-strain curve given by the Swift equation: where ε p is the effective plastic strain, and C, ε 0 and n are material parameters. The Swift curve parameters have been obtained by fitting the numerical stress-strain curve to the experimental one given in Fig. 12. Comparison of the analytical curve defined by the Swift parameters: C = 809 MPa, ε 0 = 0.0013 and n = 0.232 with the experimental stressstrain curve is shown in Fig. 14.
The numerical and experimental scaled contact forcedisplacement curves are plotted in Fig. 15. It should be noted that the numerical curve in Fig. 15 has been constructed taking h equal twice the applied displacement to the plane of the hemisphere in order to use consistent parameters with those defined in Fig. 6b. Both curves plotted in Fig. 15 are very close to each other which confirms correctness of the numer- ical model. It can be noted, however, that while the slopes of the experimental and numerical loading curves are nearly the same, the difference in the slopes of the unloading curves is significant. This will be discussed in more detail further on when parameters of the unloading model will be evaluated. Figures 16 and 17 present equivalent plastic strain and stress distributions in the contact zone at two stages of loading. It can be seen that the yielding is initiated very early, which explains why the elastic Hertzian solution has a limited validity to very low levels of loading.

Validation of the analytical contact model for loading
The analytical Storåkers contact model presented in Sec. 2 will be validated against the experimental and simulation data. Figure 18 shows the scaled force-displacement relationship according to the linear Storåkers contact model compared to the experimental and simulation curves. The linear Storåkers model is given by Eq. (15). For equal size balls made of the same perfectly plastic material, it takes the form: The initial yield stress σ Y = 214 MPa, determined in the laboratory compression test, has been taken in evaluation of the contact force according to Eq. (27). It can be seen in Fig. 18 that the plastic Storåkers model predicts the contact force which agrees quite well with the experimental and numerical results in a range of moderate deformations (in the test considered up to 5% reduction of the ball dimension along the loading direction). The loading stiffness k L for the force-displacement relationship given by Eq. (27) is expressed as follows: and it amounts to 43.26 MN/m for the given data. The relationship (27) can be transformed into the equivalent scaled form: where The decrease in stiffness manifesting itself by the deviation of the experimental and numerical curves from the Storåkers model at higher values of h/R is due to large deformations of the compressed balls.

Analysis of the unloading model
In order to have more data on the unloading, the FEM simulation has performed for a number of loading-unloadingreloading. The scaled force-displacement curve obtained in this simulation is plotted in Fig. 19. It can be seen that the linear model is quite a good approximation for unloading even for relatively small deformations. The reloading curves coincide with the unloading ones. The numerical results are consistent with the experimental results shown in Fig. 10. Similarly, the unloading slope in the numerical results gets steeper when the loading level is higher.
The unloading stiffness and its variation will be determined using the experimental and numerical results given in Figs. 10 and 19, respectively. Let us assume that the stress-type and strain-type variables, F/4R 2 and h/2R, respectively, during unloading are related linearly: where κ U is a certain modulus dependent on the achieved maximum force or deformation. Rewriting Eq. (31) in the form: and comparing it with Eq. (5) we can easily notice that the unloading stiffness k U is given by: Having calculated the unloading stiffness k U , the parameter α used in Eq. (24) can be obtained as: The   Variation of the stiffness in terms of the coefficient α is given in Fig. 21.
The plots in Figs. 20 and 21 reflect the difference between the unloading slopes in the experimental studies and numerical simulations which can be observed in Fig. 15. Experimental studies have been performed taking care to compensate properly the machine stiffness in order to eliminate its influence on the results. It was impossible, however, to avoid an error due to an indentation of the compressed balls into the plates. These problems do not occur in the FEM simulations; therefore, as far as the unloading stiffness is concerned, the numerical results are more credible than the experimental ones obtained in the tests performed in this work.
It can be further observed in Figs. 20 and 21 that the points representing the unloading moduli lie approximately along straight lines. This shows that the linear relationships postulated in Eqs. (5) and (9) are justified. The straight lines approximating the unloading moduli, however, intersect the vertical axis above the point representing the loading stiffness Analogous relationships can be written for the unloading in terms of the modulus κ unl and in terms of the scaling coefficient α Equations (35)-(37) define the respective lines plotted in Figs. 20 and 21, and their parameters k 0 U , κ 0 U and α 0 represent the initial unloading modulus, B, b and β define the variation of the unloading modulus with an increasing deformation, and k L , κ L and α L define the loading stiffness. These parameters determined from experimental and numerical data are given in Table 1.

Conclusions
The analysis of the experimental and numerical results has confirmed some known observations as well as it has given a new insight into the inelastic contact between spheres which can be useful to model the interparticle contact in the discrete element method. Both experiments and numerical results have shown that the elastic model has a very limited validity if plastic effects should be taken into account. The plastic deformation is initiated at a low load and the contact force cannot be predicted using the elastic contact models, for instance the Hertz one. Both numerical and experimental results have shown that the force-displacement curve can be approximated by a linear relationship in quite a large range of elastoplastic deformation of the contacting spheres. It has been shown that the linear Storåkers model predicts the force-displacement relationship which agrees very well with the experimental and numerical results.
Much of the attention in the present work has been paid to the modelling of the unloading in the contact between inelastic spheres. It has been demonstrated that the unloading can be approximated accurately by the linear model with the modulus varying linearly with the maximum deformation achieved in the contact. Based on the analysis of the experimental and numerical data, it has been proposed that the initial unloading stiffness can be larger than the loading stiffness. The parameters for the unloading model has been evaluated. Use of the size-independent parameters to describe the model allows us to assume that the results can also be used for other materials and particle size in the discrete element method.
The presented results and their analysis have shown that the Walton-Braun-type model with the loading stiffness evaluated according to the linear Storåkers model and the linear elastic unloading with linearly varying unloading stiffness should be an efficient and accurate model for the inelastic contact in the discrete element method using spherical particles. Obviously, one should be aware of the limitations of this model.

Conflicts of interest
The authors have declared that no conflict of interest exists.
Human and animal participants This article does not contain any studies with human participants or animals performed by any of the authors. Informed consent was obtained from all individual participants included in the study.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecomm ons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.