Contact parameter analysis and calibration in discrete element simulation of rice straw

Discrete element method was used to study and analyze the interaction between rice straws and between rice straw and agricultural machinery parts, thereby providing a scientific basis for post-harvest paddy field processing. Calibrations of rice straw-rice straw, rice straw-agricultural machinery part contact parameters (collision recovery coefficient, static friction coefficient and rolling friction coefficient) constitute an important prerequisite for the discrete element research process. In this study, the collision recovery coefficients of rice straw-steel and rice straw-rice straw were 0.230 and 0.357, respectively, which were calibrated by the collision method. The static friction coefficient and rolling friction coefficient of rice straw-steel were 0.363 and 0.208 respectively, which were calibrated by the inclined plate method and the slope method. The static friction coefficient and rolling friction coefficient of rice straw-rice straw were 0.44 and 0.07, respectively, which were calibrated by the split cylinder method. The paired t-test showed insignificant differences between calibration parameter simulation results and the physical test values (p>0.05). Taking the angle of repose that reflecting rice straw flow and friction characteristics as the evaluation index, the verification tests of the above calibration values indicated that the simulated angle of repose has no significant difference from the physical test value (p>0.05). The side plate lifting test on rice straw of different lengths showed no significant difference between the simulated angle of repose and the physical test value (p>0.05). This study can provide a basis for contact parameters choice in discrete element simulation analysis with rice straw-rice straw and rice straw-agricultural machinery parts as the research object. The calibration method can provide a reference for the contact parameter calibration of other crop straws.


Introduction 
In post-harvest rice processing, rice straw picking, bundling as well as returning are important links in modern agriculture [1][2][3] . Studying the interaction of rice straw-rice straw and rice straw-agricultural machinery parts is a basic scientific issue in the development and improvement of rice post-harvest processing equipment [4][5][6][7] . In recent years, the simulation analysis technology of discrete element method (DEM) receives great attention in the design and improvement of agricultural machinery parts. DEM can accurately analyze the movement law of agricultural materials [8,9] , provide a microscopic visual analysis of rice straw particles, dynamic meso-contact behavior between rice straw and agricultural machinery parts, thereby providing technical parameters for trial production of physical prototypes and avoiding structural problems before trial production. However, in discrete element simulation, contact parameters (collision recovery coefficient, static friction coefficient, and rolling friction coefficient) of rice straw-rice straw and rice straw-agricultural machinery parts cannot be directly determined by physical tests. The rationality of these parameters directly affects the accuracy of the simulation results [10,11] . At this stage, research on contact parameters of rice straw-rice straw and rice straw-agricultural machinery parts is still rarely reported. Domestic and foreign scholars have carried out the following research on the selection and calibration of straw contact parameters.
Fang [12] determined the wheat straw contact parameters with reference to Lenaerts et al [13] . DEM simulation was used to study the interaction mechanism of wheat straw-soil-rotary tiller. By analyzing rotary tiller force in three directions, it was concluded that increased direct straw burial in operations could reduce energy consumption in straw burial via repeated rotary tillage. Zhang [14] determine the contact parameters of rice straw-soil-rotary tiller roller by empirical values, and used DEM simulation to analyze the combined knife roller operation process in high straw returning. Liu et al. [15] selected wheat straw-steel contact parameters to establish a flexible wheat straw model using DEM and evaluate the impact of the calibrated bonding parameters on the model. The results indicated that flexible wheat straw could accurately reflect wheat straw bending property. Zhou et al. [16] cited contact parameter value in literature [14] to conduct simulation and experimental research on the spatial distribution effect after straw returning. Using the orthogonal method, Zhang et al. [17] calibrated the contact parameters between maize stalks and between maize stalks and crusher hammer, finding a relative error of 8.127% between the simulated and experimental value of the radial accumulation angle. With the angle of repose as the basis, Wu et al. [18] calibrated the discrete element parameters of tomato straw by orthogonal method to reveal the vertical spiral mixing mechanism of tomato straw.
In summary, at present, discrete element contact parameters of rice straw mostly adopt universal value or empirical value, and the large errors directly affect the accuracy of simulation analysis. To more accurately determine the contact parameters of rice straw-rice straw, rice straw-agricultural machinery parts, this study adopted virtual calibration method [19] , provided measured macro index values based on physical test, set simulation model within the measured macro index value range, and then conducted physical test and simulation test simultaneously based on specific contact parameter values of rice straw-rice straw and rice straw-agricultural machinery parts. In case of difference between the two, the target contact parameter value will be corrected and adjusted to finally determine calibrated contact parameter [20] .
That is, first, determining the contact parameter level value in the simulation test plan according to the value range of some physical test-derived contact parameters. Then, use the measured index value as the simulation test target to calculate the rice straw-steel contact parameter and rice straw-rice straw collision recovery coefficient, perform significance analysis on the simulated value and the physical test value under the calibrated values. Secondly, based on Central Composite Design (CCD), the value range of the rice straw-rice straw static friction coefficient and rolling friction coefficient derived from the pre-simulation test was used as the factor range of the simulation test. The angle of repose reflecting rice straw flow and friction characteristics is taken as the evaluation index for the simulation test. The simulation results were subject to variance analysis, and the angle of repose measured by the split cylinder method was used as the target value for optimization of the simulation regression model to get the optimized contact parameter calibration value. Finally, the verify accuracy of the above calibrated contact parameters as well as the feasibility of the calibration method; used side plate lifting test to verify universal applicability of the calibration value and the calibration method to rice straw of different lengths.

Establishment of rice straw model
The study found that the length of harvested rice straw should be no more than 100 mm [21] , and that of rice straw for other purposes should generally be no more than 30 mm [22] . The rice straw sample in this study is natural air-dried rice straw (Miaodao No. 2) taken from the Agricultural Experimental Base of Jilin University in May 2020 (125°15'E, 43°57'N). A randomized complete block design with six replicates was established in the field, and then five one-square-meter matrices were selected randomly. Considering simulation efficiency and rice straw consistency, the rice straw with an internode length of 10 mm was taken as the research object [23] . Fifty internode rice straws were randomly intercepted in each matrice (5 matrices×6 replicates, 30 in total), all rice straws were mixed evenly and sampled by the quartering method before the test. The moisture content of rice straw measured by XY-102MW halogen moisture analyzer was 10.83%. Bonding effect between rice straws can be ignored, so Hertz-Mindlin (no slip) contact model in the discrete element simulation software EDEM was selected [24] to establish the rice straw model as shown in Figure 1. The rice straw appearance observed under the ZEISS stereo microscope is shown in Figure 1a. The rice straw was approximately cylindrical with outer striped protrusions. The current straw models mostly adopt multi-balls that come with the software to form straw strips [13,14,16,18,23] , which were quite different from the actual straw in shape. In this study, the rice straw Pro/E three-dimensional model was imported into EDEM software. As shown in Figure 1b, the 10 mm long rice straw model filled with randomly generated small balls was established with model contour as the filling boundary. According to the average segment diameter of the sample, the model diameter was determined to be 5.5 mm.
a. Observation test on the outer surface of rice straw b. Single rice straw and rice straw row model Figure 1 Establishment of discrete element model of rice straw

Determination of intrinsic parameters in discrete element simulation
The material property parameters required in EDEM simulation include intrinsic parameters and contact parameters. The test process involved the contact between rice straw and rice straw, and between rice straw and collision surface. The contact material was Q235 steel commonly used in agricultural machinery [25] .
Intrinsic parameters mean characteristics of the material itself, which are usually relatively fixed. Rice straw density is measured by the liquid displacement method [26] , and intrinsic parameters required for the simulation herein are determined as listed in Table  1 through experiments and literature consultation [14,15,27,28] . Material contact parameters are related to objects in contact. The contact parameters calibrated herein include collision recovery coefficient, static friction coefficient and rolling friction coefficient of rice straw-steel and rice straw-rice straw.  [14,15] 196 a 1×10 6 [14,27] Steel 0.3 [14,15,28] 7865 [14,28] 7.9×10 10[14,28] Note: The values marked with the letter a in the table are the values measured by experiments.

Calibration test of rice straw-steel contact parameters 2.3.1 Rice straw-steel collision test
Collision recovery coefficient is a parameter that reflects the ability to return to its original shape after the collision between particles or between particles and objects. It can be expressed as the ratio of the separation speed of the two objects after the collision to the approaching speed before the collision [29,30] .
where, e is the collision recovery coefficient; 0 1 v and 0 2 v are the speeds before the collision of object 1 and object 2, respectively, m/s; v 1 and v 2 are the speeds after the collision of object 1 and object 2, respectively, m/s. Material collision recovery coefficient is commonly determined by free fall collision method or inclined plate collision method, which requires normal material rebound after falling collision [17,23,31,32] .
However, rice straw anisotropy makes accurate rebound difficult after a collision. Therefore, based on the law of conservation of kinetic energy, this study adopted a simple pendulum test similar to the Newtonian pendulum test in principle to determine the rice straw-steel collision recovery coefficient e s . The test device is shown in Figure 2. Wang [33] found that steel plate thickness has an insignificant effect on e, so the test used a steel plate of 100 mm×100 mm×8 mm. Before the test, super glue was used to fix the rice straw on one end of a nylon wire perpendicular to the rice straw axis with a diameter of 0.25 mm, while the other end was fixed on a height-adjustable support cantilever. The nylon wire has a negligible energy dissipation effect on the system [33] . The collision process photo collected by a high-speed camera (Phantom v711, Vision Research, USA) is shown in Figure 3a. The vacuum pump suction nozzle adsorbs and maintains the rice straw at the release height, while the rice straw axis is kept parallel to the horizontal plane and the steel plate plane. After turning off the vacuum pump, the rice straw was released without initial velocity to collide with the steel plate for the first time, which then rebounded to the highest point. Equation (2) of e s can be derived from Equation (1), and e s value can be calculated based on value reading on millimeter graph paper.
where, e s is the rice straw-steel collision recovery coefficient; n i v and v j are the normal rebound velocity and incident velocity of the rice straw, m/s; H a and h a are the rice straw release height and rebound height, respectively, mm; g is the acceleration of gravity, m/s 2 .
The collision simulation test is shown in Figure 3b. In the equivalent collision model, a pipeline centerline plane is a vertical plane, the pipeline section is a rectangle of 14 mm×8 mm, and the closed right end of the pipeline is the collision plane. Wang [33] found that the intrinsic parameters of the pipeline had insignificant effects on the collision process, according to the steel setting. The rice straw-steel friction coefficient was set to 0 [19] . According to EDEM user manual, the simulation time steps herein were all set to 20% of the rayleigh time step, and the grid size was 2.5 times of the minimum particle radius. The rice straw is confined in a 1/4 circular pipeline with a radius of 200 mm, and a rice straw model was generated at the vertical height H a from the lowest point of the pipeline. The rice straw was released without initial velocity, which rebounded to the highest point after colliding with the right end face of the pipeline. h a value was calculated after software post-processing.  Figure 4a. The steel plate is fixed in the groove of the self-made inclined plate instrument. To prevent rice straw rolling, the four sections of rice straw are bonded into rice straw row and placed on the steel plate. The inclined plate is slowly raised at a constant speed, the digital display goniometer displays the value of the inclination angle α 1 in real-time, and the high-speed camera captures the instant value α 1 when the rice straw row slide on the inclined plate. The tangent value is the rice straw-steel static friction coefficient, which is expressed as [23,34] : The inclined plate simulation test is shown in Figure 4b. e s was set to the calibration value in Section 2.3.1, and the rice straw-steel rolling friction coefficient was set to 0. The rice straw row shown in Figure 1d is generated on the surface of the steel plate model. The steel plate rotated counterclockwise around the Z-axis at a speed of 10°/s [19] , and the value α 1 was recorded when the rice straw row had a sliding trend.  Figure 5a. The rice straw is released from the inclined plate at a distance s 0 from the intersection line without initial velocity, which, driven by resistance, rests on the plate with distance s from the intersection. Equation (4) can be derived from the law of conservation of energy, and Equation (5) of the rice straw-steel rolling friction coefficient μ′ 1 can be derived [35,36] . In the test, to prevent rice straw bouncing due to excessively inclined angle β, or insufficient rolling due to too small β, s 0 is selected as 100 mm and β is selected as 30° through a trial test. The rolling distance s is measured and substituted into Equation (5) to get μ′ 1 .
where, μ′ 1 is the rice straw-steel rolling friction coefficient; s 0 is the distance from the initial rice straw position to the intersection line, mm; s is the rice straw rolling distance on the plate, mm; β is the slope inclination angle, (°). The slope simulation test is shown in Figure 5b. e s was set as the calibration value in Section 2.3.1, and μ 1 was set as the calibration value in Section 2.3.2. The rice straw was released where s 0 was 100 mm without initial velocity. The rice straw rests on the plate after rolling. The legend in the figure shows the real-time rolling distance of the rice straw, and the s value under static rice straw was recorded.
1. Plate 2. Intersection line 3. Inclined plate Note: s is the rice straw rolling distance on the plate, mm; s 0 is the distance from the initial rice straw position to the intersection line, mm; β is the slope inclination angle, (°).
a. Physical test b. Simulation test Figure 5 Rice straw-steel rolling friction test 2.4 Calibration test of rice straw-rice straw contact parameters 2.4.1 Rice straw-rice straw collision test Rice straw-rice straw collision recovery coefficient e r was determined by double pendulum test [33,37] with test device the same as that in Figure 2. In the test area, two nylon wires of equal length were used as the cycloid, one end of which was suspended with two rice straws approximate in shape and size. The suspension angle, height and focus point are consistent. Figure 6a shows the physical test of rice straw collision under a high-speed camera. The rice straw with initial release height reached the highest point after the first collision. Equation (6) where, e r is the rice straw-rice straw collision recovery coefficient; H 0 is the rice straw release height, mm; H b and H c are the heights of two rice straws after a collision, respectively, mm. As shown in Figure 6b, the collision simulation test adopts an equivalent collision model [34] with the basic settings the same as in section 2.3.1, in which the rice straw-steel contact parameter and the rice straw-rice straw friction coefficient are all set to 0. The rice straw is confined in the semicircular pipeline with a radius of 200 mm as shown in Figure 6b. One rice straw model is generated at the lowest point of the pipeline and the other at the distance H 0 from the vertical height. The rice straw at H 0 is released without initial velocity, and the two rice straws reach the highest points H b and H c respectively after the first collision.  The natural angle of the repose of the material reflects its flow and friction characteristics [38] . Studies have shown a significant impact on the angle of repose from friction coefficient [24,39] . In this section, by the split cylinder method, the angle of repose test was used for calibration of rice straw-rice straw static friction coefficient μ 2 and rolling friction coefficient μ′ 2 [39,40] .
The physical test is shown in Figure 7a. The steel hollow cylinder and the cylinder base have an inner diameter of 100 mm, a wall thickness of 2 mm, and a height of 65 mm. Fill the cylinder with rice straw, and slowly lift the hollow cylinder at a uniform speed. Then, determining the mutually perpendicular generatrix of the stable cone formed by the rice straw pile in four directions, and take the average value of the inclination angle as the angle of repose.
The simulation test and settings are shown in Figures 7b and  7c. The contact parameters e s , μ 1 , μ′ 1 and e r are set to the calibrated values above, μ 2 and μ′ 2 are set according to the experimental plan design. Rice straw particles are generated in the cylinder. After the particles are stable, the hollow cylinder is lifted vertically at a rate of 0.001 m/s [38] . To avoid errors caused by direct measurement of rice straw pile, as shown in Figure 8

Data analysis
In this study, Design-Expert 8.0.6.1 (Stat-Ease Inc., USA) was used for statistical analysis. A t-test was used to verify whether the simulation test results under calibrated contact parameters are significantly different from the physical test results in the collision test, inclined plate test, slope test, and side plate lifting test. F-test was used to test whether there is a statistical difference between the simulated angle of repose and the physical test value in the split cylinder method. The angle of repose in the test was read by Matlab 9.5.0 (MathWorks Inc., USA).
OriginPro 2019b (OriginLab Inc., USA) was used to plot the impact of different release heights on the collision recovery coefficient.  The measured e s values were 0.365, 0.362, 0.360, respectively. Lu et al. [29] found that H a has no significant effect on e. Therefore, the simulation test is performed by setting H a =70 mm as the initial condition, then h a =9.17 mm as the simulation target value, and the measured e s under this condition has a range of 0.335-0.392. The increments in the various test levels were determined according to the research method of Lenaerts et al. [13] . e s level value and simulation results are listed in Table 2. With rice straw-steel collision recovery coefficient as the independent variable x 1 and the rebound height as the dependent variable y a , the curve fitting equation for the simulation test is derived as: 2 11 333.14 204.12 39.576  Figure 9 shows the impact of different H a on e s . e s has a smaller simulated value than the physical test value. The error between the two was because air resistance and operation error in the physical test lead to deviation of rice straw trajectory from the simulation test, which is consistent with the reason for errors in e value of maize kernels in the study of Wang [33] . In addition, Wong et al. [37] also showed that e measurement is affected by air resistance. Rice straw outer surface has more evenly distributed strip-shaped protrusions in simulation than in physical test, which has line contact collision with the steel plate. The collision contact is larger compared to the physical test, with big energy loss, reduced rebound speed, so the simulated value is smaller than the physical test value, which is similar to the conclusion of Wang et al. [41] when analyzing e value of maize kernels in different shapes. In addition, as H a increased, e s decreased since bigger H a increased the rice straw-steel plate collision deformation. Due to extremely short collision contact time, rice straw consumes more energy in deformation and has decreased rebound speed after the collision. Therefore, e s decreased. This is similar to the reason for smaller e value under greater drop height when Yang et al. [32] calibrates e value of castor capsule and Q235 steel.  Table 3. With the rice straw-steel static friction coefficient as the independent variable x 2 and the inclination angle as the dependent variable y α1 , the curve fitting equation for the test results is The determination coefficient R 2 = 0.999. y α1 = 20.13°, so x 2 = 0.363. Set μ 1 as 0.363. According to the three sets of repeated simulation tests, α 1 was 18.14°, 17.62°, and 19.16°, respectively. Paired t-test with the physical test value showed an insignificant difference between the two (p>0.05), so the μ 1 calibration value was determined as 0.363.
The error between the μ 1 simulated value and the physical test value is due to the natural properties and uneven physical form of the rice straw in the physical test relative to the simulation model. This is similar to the study on the static friction coefficient of potatoes by Liu et al. [19] The surface roughness and viscous damping effect of rice straw in the physical test resulted in a bigger physical test value than the simulated value, which is consistent with the analysis result of DEM simulation friction coefficient of the glass bead particles calibrated by Angus et al. [42] in the shear-cell test.

Rice straw-steel rolling friction coefficient
Three sets of slope physical tests were performed, with each set repeated 5 times to take the average value. The obtained s was 149.21 mm, 151.46 mm, 150.42 mm, with an average value of 150.36 mm. The calculated μ′ 1 was 0.211. μ′ 1 level value in the simulation test can be determined based on μ′ 1 range of 0.149-0.343 in the preliminary test. The simulation results are listed in Table 4. With rice straw-steel rolling friction coefficient as the independent variable x 3 and the rolling distance as the dependent variable y s , the curve fitting equation for the test results is 2 33 3880.50 2857.80 576.52 The determination coefficient R 2 = 0.996. y s = 150.36, so x 3 = 0.208. Set μ′ 1 as 0.208, and by three sets of repeated simulation tests, resulting s is respectively 151.84 mm, 154.46 mm and 152.50 mm. Paired t-test with the physical test value showed an insignificant difference between the two (p>0.05), so μ′ 1 was calibrated as 0.208.
The error between μ′ 1 simulated value and physical test value is due to the uneven rice straw surface morphology, invisible dust to the naked eye in the physical test, which makes the physical test value greater than the simulated value. This is similar to the reason for the calibration error in DEM simulation friction coefficient of glass bead particles according to Angus et al. [42] . Steel plate surface roughness and lubrication treatment mode also explain the error between simulated value and physical test value. Ketterhagen et al. [36] also showed that the collision material surface state inevitably impacts the rolling friction coefficient. In addition, there was a minor difference between the rice straw shape and the simulation model in the physical test, which led to deviations of the rice straw rolling trajectory from the simulation. The rolling distance then becomes smaller, so μ′ 1 the physical test value was greater than the simulation value. As confirmed by Nguyen et al. [43] , the difference between soybean shape and sphere affects the grain trajectory and therefore rolling friction coefficient.  The reason for the difference between e r simulated value and The physical test value is similar to the finding of Wang [33] when measuring e r of two maize kernels. It is because the two rice straws have stable collisions in the simulation, while random error in the physical test will cause slight deflection in the rice straw collision. Rice straw in the physical test has natural attributes different from the rice straw model, resulting in errors between the two. The test found that as H 0 increased, e r decreased. The reason was that bigger H 0 increased the potential energy of the released rice straw, the deformation increased while decreasing e r . This is consistent with the conclusion of Yang et al. [32] when calibrating the collision recovery coefficient between castor capsules. 3.2.2 Rice straw-rice straw static friction coefficient and rolling friction coefficient

Analysis on calibration results of rice straw-rice straw contact parameters
The rice straw angle of repose was determined by three sets of physical tests by the split cylinder method. Each test was repeated 5 times to obtain the measured angle of repose of 34.17°, 33.35° and 32.92°, and the average value of 33.48°.
After the pre-simulation test, it was determined that μ 2 and μ′ 2 have value ranges of 0.35-0.45 and 0.06-0.08, respectively. Central Composite Design (CCD) test was performed by Design-Expert software. The test factors were the static friction coefficient A and rolling friction coefficient B between the rice straw. With test index as simulation angle of repose Y, the test factors are coded as shown in Table 6. The angle of repose is simulated with results shown in  The model's significance test has p<0.001, determination coefficient R 2 = 0.932, correction absolute coefficient R 2 adj = 0.884, which are all close to 1, indicating that the regression model can well predict the target angle of repose. With lack of fit p=0.4444> 0.05, coefficient of variation CV=0.61%, the model fits well. The regression variance analysis of the test results is shown in Table 8. The analysis shows that static friction coefficient A, quadratic term A 2 and rolling friction coefficient B have a significant effect on the angle of repose (p<0.01); quadratic term B 2 of rolling friction coefficient has a significant effect on the angle of repose (p<0.05). Ma et al. [20] confirmed that the static friction coefficient between alfalfa straws, its quadratic term and the rolling friction coefficient had a significant effect on the angle of repose, while the quadratic term of the rolling friction coefficient had a significant effect on the angle of repose. The software optimization module was used for optimization so that the simulated angle of repose was closest to the physical test value of 33.48°. Thus, μ 2 and μ′ 2 could be obtained as 0.44 and 0.07.

Analysis of contact parameter verification results
Three sets of simulation tests were performed by the split cylinder method based on contact parameter calibration values in Table 9. Each test was repeated 5 times to get average simulated angles of repose of 33.76°, 34.04°, 32.43°. Paired t-test with the physical test values showed an insignificant difference between the measured and simulated values (p>0.05). Using the side plate lifting method, Zhang et al. [17] measured radial accumulation angle of maize straw of different lengths, revealing insignificant differences between the simulation results and the physical test value. Wen et al. [24] determined that the angle of repose of urea particles with different moisture contents in the side-plate lifting test, proving the insignificant difference between the simulation and the physical test value, so it was concluded that the side-plate lifting method is a universally applicable method in the calibration of discrete element contact parameter.
Taking into account different rice straw lengths after harvest, this study adopted a side-plate lifting method to determine the angle of repose of 10-100 mm rice straw, as shown in Figure 10. As shown in Table 10, the paired t-test between the simulation test results and the physical test values showed an insignificant difference between the two (p>0.05). Also using steel side plate, Wen et al. [31] measured the angle of repose of sugarcane stalks in different lengths, confirming that the simulated value was not significantly different from the real experimental value. It indicated that the contact parameters calibrated herein have universal applicability to the discrete element simulation study of rice straw-rice straw and rice straw-agricultural machinery parts. a. Physical test b. Simulation test Figure 10 Angle of repose test by side plate lifting method

Conclusions
Based on Hertz-Mindlin (no slip) model of DEM, this study combined physical test and simulation test to calibrate rice straw-steel collision recovery coefficient, static friction coefficient and rolling friction coefficient as 0.357, 0.363 and 0.208, respectively. Rice straw-rice straw collision recovery coefficient, static friction coefficient and rolling friction coefficient were calibrated as 0.230, 0.44 and 0.07, respectively. There was an insignificant difference between the verification simulation value and physical test value under this calibration value (p>0.05). Also, no significant difference was found between the simulated and physical test value of rice straw angle of repose under different lengths (p>0.05). This study provides key contact parameters for the discrete element simulation of rice straw and its related machinery, which help with the design and optimization improvement of agricultural machinery and tools, and promotes the development process in the field of agricultural machinery.