Experimental method for identifying constitutive parameters of rolling resistance between particles in discrete element method

A method to identify the parameters of rolling resistance between particles in the discrete element method is proposed in this paper. Experiments revealed that a free rolling particle would swing back and forth marginally before it stops. We assume that the restoring force of the swing is mainly rolling resistance. An optical experimental system was designed to obtain the rolling resistance in the DEM model. A method to identify the parameters using the time history curve of the swing is proposed. We measured the time history of the angular displacement of cylindrical samples made of rubber and aluminum. In addition, the rolling stiffness and damping coefficients were identified. The identified parameters were applied to the discrete element model to simulate the experimental process. The results show that the time-history curve of swing obtained by DEM is in good agreement with the experimental curve. This verifies that the particle swings back and forth because of rolling resistance and that the parameters of the rolling resistance model can be determined by this experimental method


Introduction
The macroscopic mechanical behaviour of granular media materials is related to the microscopic characteristics of particle contact. The constitutive law for rolling resistance has a significant influence on the behaviour of granular systems simulated using the discrete element method. Bardet and Huang [1] are probably the first to introduce rotational constraints into the discrete element model. They observed that the simulation results without contact couple between the particles are outside the range of theoretical values. They consider that the contact couple may play an important role. They also demonstrated that the internal friction angle of the particle system increases when the particle rolling is fixed. Morgan [2] introduced rotational damping into a discrete element model to simulate a fault gouge. The results were close to laboratory estimates. Sakaguchi et al. [3] introduced the concept of "rolling friction" into the discrete element model and studied the plugging of granular flow in the silo discharge through experiments and numerical simulations. They observed that an arch formed by circular disks was unstable in the DEM model without rolling friction. The DEM model with rolling friction that was introduced effectively formed arches, which was observed in the physical experiment. Iwashita and Oda [4] observed large voids and high rotational gradients in a shear band experiment. However, the numerical simulation without considering rolling resistance between particles could not obtain a similar result. They determined the rolling resistance as the cause for arching. Furthermore, the rolling resistance facilitates the formation of large voids in physical experiments. To obtain numerical results close to the experimental results, they introduced rolling resistance to the conventional discrete element method and proposed a modified discrete element method (MDEM). The rolling resistance model is composed of elastic rotational springs, a rotational dash pot, a no-tensile joint, and a rotational slider. Using this model, they successfully predicted the dilatancy behaviour of the shear band, which is similar to the experimental results for natural soil. The rolling stiffness was assumed to be proportional to the contact normal force in Iwashita and Oda [4]. However, it was subsequently modified to be proportional to both contact normal force and overlap width of the contacting particles [5]. Many researchers have proposed a series of rolling resistance models inspired by Iwashita and Oda's work [6][7][8][9][10][11]. Jiang et al. [12] observed that most of the rolling resistance parameters in MDEM need to be determined by trial and error. Based on a theoretical analysis, they proposed a new definition of pure sliding and pure rolling, and introduced a "shape parameter " in the rolling resistance model. It is a dimensionless parameter that is equal to the ratio of the width of the particle contact surface to its radius. MDEM is an effective method for establishing a contact constitutive model to describe the resistance behaviour between particles.
Although the rolling resistance discrete element contact constitutive model is being continuously improved by researchers, is difficult to determine the constitutive model parameters of rolling resistance. A good experimental method to identify a significant number of parameters is unavailable. There are two reasons for this. One is that many factors influence rolling resistance, and the mechanism of rolling resistance is not clear. At present, the more widely accepted theory is that rolling resistance is affected by micro-slip on the contact surface, elastic hysteresis of the material, particle irregular shape, and contact surface adhesion [13]. The other reason is that the rolling resistance is significantly smaller than the other resistances, and it is generally accompanied by sliding resistance during the rolling process. It is difficult to experimentally determine the mechanical behaviour under pure rolling resistance. Therefore, most of the rolling resistance parameters in the DEM model are determined by trial and error.
We observed a phenomenon wherein circular particles swung back and forth in the process of free rolling. It was close to the pure rolling state, and the main restoring force was rolling resistance. An optical experimental platform with a laser displacement sensor for detecting marginal rolling of particles was designed. The time-history curves of the rolling angular displacement of cylindrical samples of various sizes were obtained using this experimental platform. The cylindrical samples were made of rubber and aluminum. A method for identifying the constitutive parameters, rotational stiffness, and rotational damping coefficients in the MDEM model is proposed based on the mechanical theory of the rolling motion of particles and the experimental results. We applied the identified parameters to the MDEM model to simulate the experimental process. The results were consistent with the experimental phenomena.

Parameters of rolling resistance in MDEM
The constitutive model of contact particles in the MDEM (Fig. 1) can be divided into three components, which are in the normal, tangential, and rolling directions. The normal contact model is composed of normal springs, dashpots, and no-tension joints. The springs are used to provide an elastic restoring force, and the dashpots are used to provide viscous resistance. The contact model in the tangential and rolling directions is similar to the model in the normal direction. All these have springs that reflect elastic behaviour and dashpots that provide viscous resistance. In addition, the model in the tangential direction has a slider that reflects the resistance caused by sliding friction. The model in the rolling direction has a rotational slider that reflects the resistance caused by the rolling friction. Jun Ai et al. [13] compared the simulation results of different rolling resistance models and observed that the MDEM model is applicable to most discrete element simulation problems.
In the MDEM, the contact rolling resistance can be described as where K r is the rolling stiffness coefficient, c r is the rolling damping coefficient, and r is the critical rolling coefficient. r expresses the critical state wherein the particles can roll forward continuously. It can be determined conveniently by physical tests. However, it is highly difficult to determine K r and c r directly by experiment. K r is generally determined from the tangential and normal stiffness by multiplying these with a proportionality coefficient [14]. Although the tangential and normal stiffness can be measured conveniently by physical experiments, the proportional coefficient is determined by trial and error. Certain experimental methods for detecting rolling resistance directly have been studied [15][16][17][18]. The rolling resistance is determined mainly based on the principle of energy conservation and the theory of static equilibrium. However, in these experimental methods, the rolling behaviour of particles is affected by rolling resistance as well as by other resistances that need to be considered. To measure the rolling resistance accurately using an experimental method, a rolling state involving only rolling resistance should be determined. We observed that a free rolling particle swung back and forth before it stopped. This is close to the pure rolling state, in which other resistances can be omitted. We suppose that the rolling resistance can be determined accurately from the swing curve of a rolling particle, and the curve can be measured using an experimental system.

Optical experimental system
A schematic of the rolling experiment of cylindrical samples is shown in Fig. 2a. An optical experimental device for detecting micro-displacements was developed. The displacement was measured by the following process: A cylindrical sample was placed on an optical experimental platform. Several contact points were selected at different locations on the circular surface of the cylindrical sample. The position of each contact point was adjusted to ensure that the sample was within the measurement range of the laser displacement sensor. We switched on the power of the experimental device, and the laser displacement sensor started to work. The laser spot position was set on the measuring point of the sample. The cylindrical sample was loaded with a marginal shock by the air impactor to swing back and forth. Meanwhile, a date acquisition terminal was used to obtain the position information of the cylindrical sample. When the entire data acquisition was completed, the displacement curve of the measured sample could be obtained by the data acquisition terminal. Next, a series of experiments were carried out in different cylindrical sample zones with various diameters, as shown in Fig. 2b and c. The displacement curves were obtained for analysis in each back and forth swing process.

Particle swing and parameter identification method
The experimental phenomenon is as follows. When circular particles cannot roll forward, these swing back and forth before stopping. The amplitude of the swing decreases gradually to zero. The back and forth swing process of particles is affected by the elastic restoring force. The swing frequency is reflected only by the rolling stiffness and rotational inertia of the particles. Because the rotational inertia can be measured conveniently, K r can be determined by the swing frequency. The attenuation of the amplitude reflects the energy dissipation during the rolling process, and c r can be determined by the attenuation of the amplitude of the swing curve. Therefore, K r and c r can be identified directly by detecting the swing displacement curve before the particle stops. Meanwhile, the influence of other resistance components on the measurement of rolling resistance parameters can be eliminated.
Because the displacement of the particle swing is highly marginal before it stops, we adopted a laser displacement sensor with high precision to measure the marginal displacement. An optical experimental device for detecting micro-displacements was established to measure the swing displacement time history data of the circular particles. The swing displacement frequency and attenuation information were obtained using an optical experimental device. Thereby, the parameters of the rolling resistance model were identified.
After impact loading, the circular particles attain a state of free rolling in the horizontal direction. The rolling process of a particle before it stops has two phases (Fig. 3a). In the first phase, it rolls forward and gradually slows down owing to the resistance. In the second phase, it cannot continue rolling forward and swings back and forth before it stops. The time-history curve of swing (Fig. 3b) could be obtained using an optical experimental device. This phenomenon was experimentally observed in this study. This back and forth swing is similar to an elastic restoring force, which can be produced only by the elastic rolling resistance. The rolling resistance of the particles can be expressed as In the back and forth swing process of the circular particle, the horizontal direction of motion is considered as the x-direction, and the rotational direction of motion is considered as the z-direction. The equilibrium equations for the x-and z-directions can be expressed as Eqs. (3) and (4): When there is no tangential deformation or relative sliding on the contact surface, the tangential damping coefficient c s = 0 , and the tangential stiffness coefficient k s = 0 . The transverse displacement x can be expressed as x = ⋅ R . The equilibrium equation for the free rolling of particles on a horizontal surface can be expressed as Eqs. (5) and (6).
If the rotational inertia of a circular particle to the center is J , the equivalent rotational inertia I g of the particle in the state of the horizontal back and forth swing can be expressed as Eq. (7). Thus, the equilibrium in Eq. (6) can be expressed using Eq. (8). In addition, the solution of Eq.(8) can be expressed as Eq. (9): The free vibration curve of the low damping system (Eq. (9)) was used to fit the swing curve of the cylindrical sample. According to Eq. (10), K r can be obtained by measuring the period T and frequency of the curve. c r can be obtained using Eq. (11) and the attenuation change in the curve:

Experimental data processing
For the swing curve that is obtained (Fig. 3b), the position of the cylindrical sample when it is completely stationary is set as the origin. The angular displacement curve is transformed using the transverse displacement curve according to the radius of the sample. The adjusted swing curve is shown by the blue curve in Fig. 4.
The free vibration curve of the low-damping system was used to fit the swing curve of the cylindrical sample. The fitted curve (fitted with Eq. (12)) is indicated by the red dashed line in Fig. 4. This figure shows that the red fitted curve is in good agreement with the blue swing curve. Therefore, the swing curve of the cylindrical sample can be processed as a low-damping vibration curve: The rotational inertia I g of the cylindrical sample can be obtained using Eq. (7). The equivalent rolling stiffness of the cylindrical sample in the swing process can be obtained using Eq. (10) as The damping coefficient is: The parameters K r and c r at different positions were obtained by processing the swing curves at other positions. A series of cylindrical samples with different diameters and materials was tested using the designed optical experimental system. The rolling resistance parameters were identified using this method.
= e −0.479t × 0.174sin(7.565t + 0.716) Figures 5 and 6 show that the rolling resistance parameters of the cylindrical samples are related to their geometric size and elastic modulus E of the material. The elastic modulus of rubber samples E rub is 6.1 × 10 6 Pa and that of the aluminum sample E al is 6.9 × 10 10 Pa. For cylindrical sample of identical material, K r and c r increase with an increase in the diameter d. As shown in Fig. 7, where K r is quadratic relation of d, which is consistent with the expression in reference [13]. c r increases linearly with d. Similarly, for cylindrical samples of equal diameter, K r and c r increase with the increase in E. According to Eqs. (10) and (11) and experimental data, . Therefore, K r and c r increase nonlinearity with E.
Poisson's ratio has little effect on the measured rolling resistance parameters. The identified parameters K r and c r display a certain degree of dispersion caused by measurement errors and uncertain factors in the experiment. Their average values are expressed as K r and c r , the standard deviation are S Kr and S cr .

Verification of rolling stiffness coefficient K r
When the angular displacement of the cylindrical sample attains its maximum, the angular velocity is zero. The point of the maximum angular displacement is set as the feature point of the test curve, as shown in Fig. 8. According to Eq. (2), the rolling resistance of the four feature points can be expressed by Eq. (15): The resistance moment of a cylindrical sample in motion can be expressed by Eq. (16): where I g represents the rotational inertia and represents the angular acceleration during the swing process. The rolling stiffness parameter K r that is identified can be verified by comparing the differences between the two calculation methods in Eqs. (15) and (16).
As shown in Fig. 8, the feature point A in the test curve is considered as an example. The curve near the peak point is polynomially fitted with Eq. (17). The second-order derivative in Eq. (18) can be obtained by fitting Eq. (17).
The angular acceleration for the feature point A is calculated according to the horizontal coordinate and Eq. (18). The rolling resistance moment is calculated according to (17) = 0.256t 3 − 4.042t 2 + 20.493t − 33.714 Eqs. (15) and (16), respectively. K r is verified by comparing the values of the resistance moment obtained by the two methods. The verification results of K r are presented in Table 1. Table 1 shows that there is a certain difference between the resistance moment calculated using Eq. (15) ( M 1 ) and and that calculated using Eq. (16) ( M 2 ). The average relative  Fig. 3 Swing state of a circular particle before it stops rolling: a swing in rolling process; b swing transverse displacement curve difference is calculated as =11.7%. Thus, the relative difference between the two methods is marginal and satisfactory. It is feasible to calculate the rolling resistance using K r , which was identified by the dynamic experiment. In conclusion, the parameter K r identified by the dynamic test can be used to establish a rolling resistance model.

Verification of rolling resistance model parameters by discrete element simulation
The K r and c r in the rolling process of the circular particles were identified by the dynamic test. Then, we substituted the identified parameters into the DEM with Eq. (1). As  Fig. 9a, the rolling process of a circular particle is simulated using the discrete element method. In this simulation, the diameter d of the particle is 80 mm, and its elastic modulus E is 6.1 × 10 6 Pa . The calculated swing curve is shown as the green curve in Fig. 9b. The rolling resistance parameters identified by the dynamic test method were substituted into the discrete element model. Fig. 9b shows that the particle simulation motion curve is in good agreement with the test curve. Therefore, the rolling resistance parameters identified by the method proposed in this paper are applicable to the DEM simulation of particle material motion.

Summary and conclusions
To establish a discrete element model of particle rolling resistance, a high-precision experimental device for measuring micro-displacement was developed, and a method to identify the rolling resistance parameters of particles (1) We observed that a free rolling particle would swing back and forth before it stops. This is because of the elastic rolling resistance. The parameters of the rolling resistance model in the DEM can be determined by this experimental phenomenon.
(2) An optical experimental platform was developed that can be used to measure the time-history curve of particle swing. In addition, we propose a method to identify the rolling stiffness K r and rolling damping parameter c r based on the time-history curve.
(3) The results of the DEM simulation were in good agreement with the experimental results. The rolling Fig. 8 Verification of rolling stiffness coefficient K r at feature points of angular displacement curve Fig. 9 Verification of rolling process of a single circular particle by discrete element simulation: a simulation of single particle rolling process by discrete element method; b comparison between the calcu-lation results of discrete element and the test curve. (parameters identified K r = 2.492 × 10 −2 (N ⋅ m/rad) , c r = 4.153 × 10 −4 (N ⋅ m ⋅ s/rad)) resistance parameters identified by the method proposed in this study are applicable to DEM simulations.