Soil Particle Modeling and Parameter Calibration Based on Discrete Element Method

: In order to establish a Discrete Element Method (DEM) model of soil particles, the soil in the laboratory soil bin was used as the research object. The soil texture was determined to be sandy loam by sieving, and the shape of the soil particles was analyzed by an image particle analyzer to establish a geometric model of the soil particles. The Edinburgh Elasto-Plastic Adhesion (EEPA) model was chosen as the contact model for the soil particle simulation analysis, and the accuracy of the model selection was determined by texture tests. The parameters in the contact model played a crucial role in the results of the simulation. Test methods were used to obtain parameters for the soil particles that were easy to measure. For parameters that could not be measured in the contact model, a direct shear test was used as the calibration test, and after screening the sensitive parameters using the PB test, the response surface method was used to calibrate the sensitive parameters. The accuracy of the calibration results was veriﬁed by comparing the simulation and test results of the direct shear test under different loadings.


Introduction
The DEM is now widely used in the field of agricultural engineering [1][2][3].During sowing, field management, and harvesting, there is inevitably interaction between soil particles, between soil particles and seeds, and between soil particles and boundaries [4][5][6].When the DEM is used to analyze the above-mentioned interactions, a DEM model of the soil particles is first established.
Ucgul et al. established a spherical model of soil particles with a radius of 10 mm and analyzed the soil after plough tillage.The results showed that the simulated values of the lateral and forward movement of the surface soil were larger than the measured values, and the simulation of the deep soil agreed well with the test results [7].Bravo et al. developed a spherical model of soil particles to simulate a soil direct shear test.Particle models with different particle sizes (3-4 mm in the surface layer, 4-6 mm in the middle layer, and 6-8 mm in the bottom layer) were used to fill the soil bin, and soil tillage in hard-dry, soft-wet, and friable states was simulated and analyzed [8].Our analysis showed that the particles in the soil particle models established by various scholars are large and spherical, which is different from the actual size and shape of soil particles.Therefore, the size and shape of soil particles need to be studied in depth.
The different textures of soil particles require different contact models for simulation analysis.Chen et al. conducted a deep loosening test and simulation on three different soils (coarse sand, loam, and sandy loam) using a parallel bond model (PBM).The simulation and test results for the soil cutting force and soil disturbance characteristics caused by deep loosening were compared and analyzed.The results showed that the relative error was less than 10% in most cases [9].Xiang et al. used the Hertz-Mindlin with JKR model (JKR model) to analyze clay loam soils in southern China.The parameters of the JKR model were calibrated by angle of repose tests.The accuracy of the DEM model for reflecting the physical and mechanical properties of the soil was verified by comparing the simulation of hole formation with the soil bin test [10].Our analysis showed that there are a variety of contact models applicable to the simulation of soil particles.However, most scholars do not provide an explanation of why the corresponding contact model was used for the test soil.Therefore, the contact model chosen for the test soil and the accuracy of the model require further analysis and research.
Soil particles have a number of parameters which can play a very significant role in simulation results.Ucgul et al. determined the DEM parameters required to simulate soil tillage processes by means of an angle of repose test and a cone penetration test.The results showed that when the appropriate contact model and parameters were adopted, the DEM could effectively predict the interaction force between non-cohesive soil and tools and simulate the movement of soil particles [7].Bravo et al. obtained the macroscopic parameter values of soil mechanical properties by a direct shear test.The values of parameters such as the elastic modulus, shear strength, friction coefficient, and cohesion of the soil were obtained by statistical regression equations.The accuracy of the model was also verified by comparing the simulation with test results using a soil bin tillage test [8].The combinations of parameters measured by the abovementioned scholars without analyzing their sensitivity are not accurate enough and may differ significantly from the actual values.When there are many parameters to be determined, the parameters that need to be calibrated and the calibration methods should be studied in depth.
Based on the abovementioned problems, a geometric model of soil particles was established in this paper after conducting a shape analysis.According to the characteristics of the test soil, the EEPA model was chosen for the contact model.The direct shear test was used as a calibration test, and sensitive parameters were screened by PB tests for unmeasurable simulation parameters.Further the sensitive parameters were calibrated using the response surface method, and the accuracy of the calibrated parameters was verified.

Soil Particle Size Distribution and Texture
The test soil was sieved through a 3 mm diameter soil sieve.A Chinese standard soil sieve was selected, as shown in Figure 1.The soil was sieved and tested with sieve diameters of 2 mm, 1 mm, 0.5 mm, 0.25 mm, and 0.075 mm, arranged in order from top to bottom, with a collection box at the base.soils (coarse sand, loam, and sandy loam) using a parallel bond model (PBM).The simulation and test results for the soil cutting force and soil disturbance characteristics caused by deep loosening were compared and analyzed.The results showed that the relative error was less than 10% in most cases [9].Xiang et al. used the Hertz-Mindlin with JKR model (JKR model) to analyze clay loam soils in southern China.The parameters of the JKR model were calibrated by angle of repose tests.The accuracy of the DEM model for reflecting the physical and mechanical properties of the soil was verified by comparing the simulation of hole formation with the soil bin test [10].Our analysis showed that there are a variety of contact models applicable to the simulation of soil particles.However, most scholars do not provide an explanation of why the corresponding contact model was used for the test soil.Therefore, the contact model chosen for the test soil and the accuracy of the model require further analysis and research.
Soil particles have a number of parameters which can play a very significant role in simulation results.Ucgul et al. determined the DEM parameters required to simulate soil tillage processes by means of an angle of repose test and a cone penetration test.The results showed that when the appropriate contact model and parameters were adopted, the DEM could effectively predict the interaction force between non-cohesive soil and tools and simulate the movement of soil particles [7].Bravo et al. obtained the macroscopic parameter values of soil mechanical properties by a direct shear test.The values of parameters such as the elastic modulus, shear strength, friction coefficient, and cohesion of the soil were obtained by statistical regression equations.The accuracy of the model was also verified by comparing the simulation with test results using a soil bin tillage test [8].The combinations of parameters measured by the abovementioned scholars without analyzing their sensitivity are not accurate enough and may differ significantly from the actual values.When there are many parameters to be determined, the parameters that need to be calibrated and the calibration methods should be studied in depth.
Based on the abovementioned problems, a geometric model of soil particles was established in this paper after conducting a shape analysis.According to the characteristics of the test soil, the EEPA model was chosen for the contact model.The direct shear test was used as a calibration test, and sensitive parameters were screened by PB tests for unmeasurable simulation parameters.Further the sensitive parameters were calibrated using the response surface method, and the accuracy of the calibrated parameters was verified.

Soil Particle Size Distribution and Texture
The test soil was sieved through a 3 mm diameter soil sieve.A Chinese standard soil sieve was selected, as shown in Figure 1.The soil was sieved and tested with sieve diameters of 2 mm, 1 mm, 0.5 mm, 0.25 mm, and 0.075 mm, arranged in order from top to bottom, with a collection box at the base.After sieving, the soil particles were weighed in the different particle size ranges, and three replicate tests were carried out to analyze the test data, as shown in    The mass percentage of soil particles was 79.39%, 9.55%, and 11.06% for the particle size ranges of 0~1 mm, 1~2 mm, and 2~3 mm, respectively.By comparing the soil texture classification criteria [11], it was clear that the soil type selected for testing was sandy loam.

Geometric Model of Soil Particles
In order to study the shape of the soil particles and facilitate the establishment of their geometric model, the soil particles were scanned and analyzed using an XS-2100 image particle analyzer, as shown in Figure 2.
After sieving, the soil particles were weighed in the different particle size ranges, and three replicate tests were carried out to analyze the test data, as shown in Table 1.The mass percentage of soil particles was 79.39%, 9.55%, and 11.06% for the particle size ranges of 0~1 mm, 1~2 mm, and 2~3 mm, respectively.By comparing the soil texture classification criteria [11], it was clear that the soil type selected for testing was sandy loam.

Geometric Model of Soil Particles
In order to study the shape of the soil particles and facilitate the establishment of their geometric model, the soil particles were scanned and analyzed using an XS-2100 image particle analyzer, as shown in Figure 2. The test procedure was as follows: soil particles with a particle size range of 0~0.075 mm were placed on a slide of the image particle analyzer.A sample of 30 soil particles was taken, and the shape of each particle was observed on a computer after magnification.Using the same method, enlarged images of soil particles with particle sizes ranging from 0.075 to 0.25 mm, 0.25 to 0.5 mm, and 0.5 to 1 mm were obtained.Figure 3a-d show twodimensional scanned images of representative soil particles of different particle size ranges selected for this paper.For soil particles in the particle size range of 1~2 mm and 2~3 mm, photographs were taken using a camera with 5× magnification, as shown in Figure 3e,f.The analysis in Figure 3 shows that the soil particles could be approximated as triangle-like and sphere-like in the different particle size ranges, and that both shapes accounted for the same proportion of the total mass.The test procedure was as follows: soil particles with a particle size range of 0~0.075 mm were placed on a slide of the image particle analyzer.A sample of 30 soil particles was taken, and the shape of each particle was observed on a computer after magnification.Using the same method, enlarged images of soil particles with particle sizes ranging from 0.075 to 0.25 mm, 0.25 to 0.5 mm, and 0.5 to 1 mm were obtained.Figure 3a-d show two-dimensional scanned images of representative soil particles of different particle size ranges selected for this paper.For soil particles in the particle size range of 1~2 mm and 2~3 mm, photographs were taken using a camera with 5× magnification, as shown in Figure 3e,f.The analysis in Figure 3 shows that the soil particles could be approximated as triangle-like and sphere-like in the different particle size ranges, and that both shapes accounted for the same proportion of the total mass.
Geometric models of the two soil particle shapes described above were established.A large number of soil particles were required for simulation, so the number of constituent spheres in the particle model had to be reduced as much as possible.
For the triangle-like soil particles, three constituent spheres of the same radius (r) were used to model the geometry.The maximum distance H between the boundaries of the constituent spheres was defined as the current particle size R of the soil particles (H = R = 3r), as shown in Figure 4. Geometric models of the two soil particle shapes described above were established.A large number of soil particles were required for simulation, so the number of constituent spheres in the particle model had to be reduced as much as possible.
For the triangle-like soil particles, three constituent spheres of the same radius (r) were used to model the geometry.The maximum distance H between the boundaries of the constituent spheres was defined as the current particle size R of the soil particles (H = R = 3r), as shown in Figure 4.For the sphere-like soil particles, the geometric model was built directly using a single sphere with the same radius as the soil particles (R = r), as shown in Figure 5.  Geometric models of the two soil particle shapes described above were established.A large number of soil particles were required for simulation, so the number of constituent spheres in the particle model had to be reduced as much as possible.
For the triangle-like soil particles, three constituent spheres of the same radius (r) were used to model the geometry.The maximum distance H between the boundaries of the constituent spheres was defined as the current particle size R of the soil particles (H = R = 3r), as shown in Figure 4.For the sphere-like soil particles, the geometric model was built directly using a single sphere with the same radius as the soil particles (R = r), as shown in Figure 5.For the sphere-like soil particles, the geometric model was built directly using a single sphere with the same radius as the soil particles (R = r), as shown in Figure 5.

Particle Population Modeling
Following the particle size distribution of soil particles in the previous section, the mass ratio of soil particles with particle sizes ranging from 0 to 1 mm, 1 to 2 mm, and 2 to 3 mm could be approximated as 4:1:1.Therefore, when modeling the soil particle population, soil particles of different particle size ranges were generated in the proportions described above.A population of soil particles in the same particle size range was generated according to a uniform distribution using the middle size of each interval as the particle radius.In this way, the distribution of the modeled soil particle population was closer to that of the actual soil particle population.
The soil particles were small in size; however, the number of soil particles required for the simulation was relatively large.In order to make the soil particle model resemble actual soil more closely, and taking into account the duration and cost of the simulation, the particle diameter was chosen to be 1 mm, and a 1-sphere and 3-sphere model were generated in equal proportions so that the parameters of the simulation as well as the simulation results would be closer to the test values.

Measurement of Soil Moisture Content
The moisture content of the test soil was within 25% ± 1%, and, in accordance with the above test requirements, soil with a corresponding moisture content was configured in this study.
The moisture content measured by the XY-102MW halogen moisture analyzer is shown in Figure 6.The moisture content of the soil sample was measured by taking 3-5 g of the prepared soil and placing it on the tray of the moisture analyzer to be heated and dried.Three replicate tests were carried out, and the moisture content of the prepared soil was found to be 24.32%, which met the test requirements.

Measurement of Soil Particle Density
The type of density to be obtained for the EDEM simulation was the density of the soil particles, which was measured using the pycnometer method [12][13][14][15].The test device is shown in Figure 6.The test was repeated three times to obtain a soil particle density of 1.841 g/cm 3 .

Particle Population Modeling
Following the particle size distribution of soil particles in the previous section, the mass ratio of soil particles with particle sizes ranging from 0 to 1 mm, 1 to 2 mm, and 2 to 3 mm could be approximated as 4:1:1.Therefore, when modeling the soil particle population, soil particles of different particle size ranges were generated in the proportions described above.A population of soil particles in the same particle size range was generated according to a uniform distribution using the middle size of each interval as the particle radius.In this way, the distribution of the modeled soil particle population was closer to that of the actual soil particle population.
The soil particles were small in size; however, the number of soil particles required for the simulation was relatively large.In order to make the soil particle model resemble actual soil more closely, and taking into account the duration and cost of the simulation, the particle diameter was chosen to be 1 mm, and a 1-sphere and 3-sphere model were generated in equal proportions so that the parameters of the simulation as well as the simulation results would be closer to the test values.

Measurement of Soil Moisture Content
The moisture content of the test soil was within 25% ± 1%, and, in accordance with the above test requirements, soil with a corresponding moisture content was configured in this study.
The moisture content measured by the XY-102MW halogen moisture analyzer is shown in Figure 6.The moisture content of the soil sample was measured by taking 3-5 g of the prepared soil and placing it on the tray of the moisture analyzer to be heated and dried.Three replicate tests were carried out, and the moisture content of the prepared soil was found to be 24.32%, which met the test requirements.

Direct Shear Test
In this paper, the geometric model of soil particles was verified and the simulation parameters were calibrated by a direct shear test [16,17].A ZJ-1B strain-controlled direct shear apparatus, as shown in Figure 7, which is an automatically controlled direct shear

Measurement of Soil Particle Density
The type of density to be obtained for the EDEM simulation was the density of the soil particles, which was measured using the pycnometer method [12][13][14][15].The test device is shown in Figure 6.The test was repeated three times to obtain a soil particle density of 1.841 g/cm 3 .

Direct Shear Test
In this paper, the geometric model of soil particles was verified and the simulation parameters were calibrated by a direct shear test [16,17].A ZJ-1B strain-controlled direct shear apparatus, as shown in Figure 7, which is an automatically controlled direct shear apparatus, was used to carry out the direct shear test of the soil.

Direct Shear Test
In this paper, the geometric model of soil particles was verified and the simulation parameters were calibrated by a direct shear test [16,17].A ZJ-1B strain-controlled direct shear apparatus, as shown in Figure 7, which is an automatically controlled direct shear apparatus, was used to carry out the direct shear test of the soil.The shear strength (τ) of the soil under different vertical loadings was measured by a direct shear test.The cohesion (c) and the internal friction angle (φ) were calculated as follows: The test procedure was as follows.We applied a vertical loading of 100 kPa and a shear speed of 2.4 mm/min and recorded the percent meter data every 30 s.When the percent meter was stable, we closed the forward button and recorded the data.The maximum shear strength ( max τ ) was calculated by the following formula: where C is the coefficient of the force ring, Rf is the reading of the percent meter when it is stable, and A is the cross-sectional area of the direct shear box.Three soil specimens were taken under vertical loadings of 100 kPa, 200 kPa, 300 kPa, and 400 kPa, and three tests were conducted for each.An image of the variation in shear strength with displacement for different vertical loading cases is shown in Figure 8.The shear strength (τ) of the soil under different vertical loadings was measured by a direct shear test.The cohesion (c) and the internal friction angle (ϕ) were calculated as follows: The test procedure was as follows.We applied a vertical loading of 100 kPa and a shear speed of 2.4 mm/min and recorded the percent meter data every 30 s.When the percent meter was stable, we closed the forward button and recorded the data.The maximum shear strength (τ max ) was calculated by the following formula: where C is the coefficient of the force ring, R f is the reading of the percent meter when it is stable, and A is the cross-sectional area of the direct shear box.Three soil specimens were taken under vertical loadings of 100 kPa, 200 kPa, 300 kPa, and 400 kPa, and three tests were conducted for each.An image of the variation in shear strength with displacement for different vertical loading cases is shown in Figure 8.
The relationship between the shear strength of the soil and the vertical loading is shown in Figure 9.The functional relationship can be expressed as: The intersection of the straight line with the vertical coordinate represents the cohesive force (c) of the soil particles, the angle between the straight line and the horizontal coordinate is the angle (ϕ) of the internal friction of the soil, and the tangent of this angle is the slope of the corresponding straight line, as shown in Figure 9. Based on Equation (3), the cohesion and internal friction angle of the test soil were calculated to be 5.94 kPa and 5.43 • , respectively.The relationship between the shear strength of the soil and the vertical loading is shown in Figure 9.The functional relationship can be expressed as: The intersection of the straight line with the vertical coordinate represents the cohesive force (c) of the soil particles, the angle between the straight line and the horizontal coordinate is the angle (φ) of the internal friction of the soil, and the tangent of this angle is the slope of the corresponding straight line, as shown in Figure 9. Based on Equation (3), the cohesion and internal friction angle of the test soil were calculated to be 5.94 kPa and 5.43°, respectively.

Simulation Analysis
EDEM software was used to simulate the direct shear test.Figure 10 shows the geometric model of the DEM simulation.At the same time, the soil particles were modeled as a 1-sphere model and a 3-sphere model with a particle size of 1 mm.After generating soil particles in the straight shear box, simulation analysis was performed, and the steps are shown in detail below.The relationship between the shear strength of the soil and the vertical loading is shown in Figure 9.The functional relationship can be expressed as: The intersection of the straight line with the vertical coordinate represents the cohesive force (c) of the soil particles, the angle between the straight line and the horizontal coordinate is the angle (φ) of the internal friction of the soil, and the tangent of this angle is the slope of the corresponding straight line, as shown in Figure 9. Based on Equation (3), the cohesion and internal friction angle of the test soil were calculated to be 5.94 kPa and 5.43°, respectively.

Simulation Analysis
EDEM software was used to simulate the direct shear test.Figure 10 shows the geometric model of the DEM simulation.At the same time, the soil particles were modeled as a 1-sphere model and a 3-sphere model with a particle size of 1 mm.After generating soil particles in the straight shear box, simulation analysis was performed, and the steps are shown in detail below.

Simulation Analysis
EDEM software was used to simulate the direct shear test.Figure 10 shows the geometric model of the DEM simulation.At the same time, the soil particles were modeled as a 1-sphere model and a 3-sphere model with a particle size of 1 mm.After generating soil particles in the straight shear box, simulation analysis was performed, and the steps are shown in detail below.

Contact Model for Soil Particles
The four main contact models involved in soil particle simulation calculations are as follows.The Hertz-Mindlin model is suitable for non-compressible and non-sticky soil and materials such as dry gravel and sand.The Hertz-Mindlin with JKR model is suitable for non-compressible, sticky, or very sticky soil and materials such as wet gravel.The hysteretic spring model is suitable for compressible or very compressible dry soil and materials such as freshly ploughed soil.The EEPA model is suitable for compressible sticky or very sticky soil or soft, very sticky soil and materials such as clay and very wet sand [18][19][20][21][22].
The test soil was a sandy loam with a moisture content of 25% ± 1%.The properties of the selected test soils were analyzed using a texture analyzer, as shown in Figure 12.The test procedure was as follows.The texture analyzer was calibrated, the test soil was placed in an aluminum foil soil tray, and the probe was moved downwards after starting the texture analyzer.After the probe had made contact with the soil, the force on the probe gradually increased.When the force reached its maximum, the moving probe started to travel in the opposite direction until it stopped.

Contact Model for Soil Particles
The four main contact models involved in soil particle simulation calculations are as follows.The Hertz-Mindlin model is suitable for non-compressible and non-sticky soil and materials such as dry gravel and sand.The Hertz-Mindlin with JKR model is suitable for non-compressible, sticky, or very sticky soil and materials such as wet gravel.The hysteretic spring model is suitable for compressible or very compressible dry soil and materials such as freshly ploughed soil.The EEPA model is suitable for compressible sticky or very sticky soil or soft, very sticky soil and materials such as clay and very wet sand [18][19][20][21][22].
The test soil was a sandy loam with a moisture content of 25% ± 1%.The properties of the selected test soils were analyzed using a texture analyzer, as shown in Figure 12.The test procedure was as follows.The texture analyzer was calibrated, the test soil was placed in an aluminum foil soil tray, and the probe was moved downwards after starting the texture analyzer.After the probe had made contact with the soil, the force on the probe gradually increased.When the force reached its maximum, the moving probe started to travel in the opposite direction until it stopped.
The relationship between the force and the displacement of the probe is shown in Figure 13.Our analysis showed that when the displacement increased in the forward direction, the force gradually increased until it reached its maximum value.Conversely, when the displacement decreased, the force gradually decreased to zero, followed by the force increasing to its maximum and then gradually decreasing until it reached a relatively stable state.The above test analysis demonstrated that the test soil was a compressible and sticky material.The EEPA model was therefore used as the contact model for the test soil simulation in order to simulate the movement of the soil particles.
Figure 13.Our analysis showed that when the displacement increased in the forward direction, the force gradually increased until it reached its maximum value.Conversely, when the displacement decreased, the force gradually decreased to zero, followed by the force increasing to its maximum and then gradually decreasing until it reached a relatively stable state.The above test analysis demonstrated that the test soil was a compressible and sticky material.The EEPA model was therefore used as the contact model for the test soil simulation in order to simulate the movement of the soil particles.

Introduction to the Simulation Parameters
The three categories of simulation parameters were as follows: 1.The material parameters included the constituent spherical coordinates, density (ρ), shear modulus (G), and Poisson's ratio (ν) of the soil particles and the density (ρ0), shear modulus (G0), and Poisson's ratio (λ0) of the boundary material.2. The parameters of interaction included the coefficient of static friction (μ1), coefficient of rolling friction (μ2), and coefficient of restitution (e) between soil particles, and the coefficient of static friction (μ11), coefficient of rolling friction (μ21), and coefficient of restitution (e1) between soil particles and boundaries.3. The parameters of the EEPA model included the contact pull-off force (f0), surface energy (γ), contact plastic ratio (λ), slope exp (n0), tensile exp (n), and tangential stiff multiplier (k). Figure 13.Our analysis showed that when the displacement increased in the forward direction, the force gradually increased until it reached its maximum value.Conversely, when the displacement decreased, the force gradually decreased to zero, followed by the force increasing to its maximum and then gradually decreasing until it reached a relatively stable state.The above test analysis demonstrated that the test soil was a compressible and sticky material.The EEPA model was therefore used as the contact model for the test soil simulation in order to simulate the movement of the soil particles.Figure 13.The relationship between force and displacement of the probe.

Introduction to the Simulation Parameters
The three categories of simulation parameters were as follows: 1.The material parameters included the constituent spherical coordinates, density (ρ), shear modulus (G), and Poisson's ratio (ν) of the soil particles and the density (ρ0), shear modulus (G0), and Poisson's ratio (λ0) of the boundary material.2. The parameters of interaction included the coefficient of static friction (μ1), coefficient of rolling friction (μ2), and coefficient of restitution (e) between soil particles, and the coefficient of static friction (μ11), coefficient of rolling friction (μ21), and coefficient of restitution (e1) between soil particles and boundaries.3. The parameters of the EEPA model included the contact pull-off force (f0), surface energy (γ), contact plastic ratio (λ), slope exp (n0), tensile exp (n), and tangential stiff multiplier (k).
Figure 13.The relationship between force and displacement of the probe.

Calibration and Verification of Simulation Parameters 6.1. Introduction to the Simulation Parameters
The three categories of simulation parameters were as follows: 1.

2.
The parameters of interaction included the coefficient of static friction (µ 1 ), coefficient of rolling friction (µ 2 ), and coefficient of restitution (e) between soil particles, and the coefficient of static friction (µ 11 ), coefficient of rolling friction (µ 21 ), and coefficient of restitution (e 1 ) between soil particles and boundaries.
The rest of the parameters (µ 1, µ 2, e, γ, λ, n, and k) were calibrated by comparing the direct shear test with the simulation, and the response index was the shear strength.

Parameter Screening
The Plackett-Burman design (PBD) test was used to analyze the sensitivity of the parameters and to select the parameters that were sensitive to the response index.The PBD test is shown in Table 2.The corresponding parameters were input into the EDEM software for simulation.Each set of parameters was simulated five times.The factor analysis was carried out using the shear strength (τ) as the response index.Based on the analysis of variance of the PBD test, it could be seen that the p-value for the coefficient of static friction was 0.031, indicating that the coefficient of static friction was significant to the response index.The p-value for the surface energy was 0.184, and the p-value for the coefficient of rolling friction was 0.14, as shown in Table 3.These two parameters were approximately significant to the response index.The p-values for the other parameters were too large to have a significant effect on the response index.The Pareto diagram of the standardization effect of the PB test was analyzed, as shown in Figure 14.The histograms above the red dashed line indicate a significant effect on the response index.It is clear from the graph that the coefficient of static friction had a significant effect on the response index.The results for the coefficient of rolling friction and surface energy were closer to the red dashed line and had an approximately significant effect on the response index.
Based on the above analysis, the parameters to be calibrated in this paper were identified as the coefficient of static friction, the coefficient of rolling friction, and the surface energy.
on the response index.It is clear from the graph that the coefficient of static friction had a significant effect on the response index.The results for the coefficient of rolling friction and surface energy were closer to the red dashed line and had an approximately significant effect on the response index.
Based on the above analysis, the parameters to be calibrated in this paper were identified as the coefficient of static friction, the coefficient of rolling friction, and the surface energy.

Parameter Calibration
The response surface method was used to calibrate the above three parameters, and the other parameters were taken as follows: contact plastic ratio, λ = 0.35; coefficient of restitution, e = 0.5; tensile exp, n = 1; and tangential stiff multiplier, k = 1.
A three-factor Box-Behnken design (BBD) test was developed, as shown in Table 4.The response index was the shear strength, with a simulated vertical loading of 100 kPa.

Parameter Calibration
The response surface method was used to calibrate the above three parameters, and the other parameters were taken as follows: contact plastic ratio, λ = 0.35; coefficient of restitution, e = 0.5; tensile exp, n = 1; and tangential stiff multiplier, k = 1.
A three-factor Box-Behnken design (BBD) test was developed, as shown in Table 4.The response index was the shear strength, with a simulated vertical loading of 100 kPa.The corresponding parameters were entered into EDEM software for simulation.Each set of parameters was simulated five times, and the simulation results were analyzed.The analysis of variance for the BBD test showed that, for the model as a whole, with a p-value of 0.001, the predictive model of the test was extremely significant.For the linear analysis, the p-value of 0 for the coefficient of static friction had a significant effect on the response index.For the squared analysis, it could be seen that µ 1 × µ 1 had a p-value of 0.001, which had a significant effect on the response index.For the two-factor interaction analysis, it could be seen that µ 1 × µ 2 had a p-value of 0.031, which had a significant effect on the response index as shown in Table 5.The analysis of the model for the BBD test showed that the R-sq value reached 98.44%, indicating a very good fit, as shown in Table 6.The corresponding regression equation was further obtained, and the formula was as follows: In the direct shear test, when the vertical loading was 100 kPa, the corresponding shear strength was 16.25 kPa.Firstly, the range of parameters was determined, wherein the coefficient of static friction was 0.1~0.9, the coefficient of rolling friction was 0.1~0.7,and the surface energy was 0~1 J/m 2 .Next, the parameters were optimized for the constrained range, and the final values of 0.9, 0.7, and 1 J/m 2 were calibrated for the coefficient of static friction, the coefficient of rolling friction, and the surface energy, respectively.

Validation of Calibration Parameters
To verify the accuracy of the calibration parameters, the shear strength of the soil was simulated and calculated for vertical loadings of 100 kPa, 200 kPa, 300 kPa, and 400 kPa.The relationship between the shear strength and the simulation time of the soil particles for the direct shear simulation is shown in Figure 15.When the shear box started to move, shear strength was generated, which gradually increased and started to decrease after reaching a maximum value, and the simulation recorded the shear strength under different loadings.
To verify the accuracy of the calibration parameters, the shear strength of the soil was simulated and calculated for vertical loadings of 100 kPa, 200 kPa, 300 kPa, and 400 kPa.The relationship between the shear strength and the simulation time of the soil particles for the direct shear simulation is shown in Figure 15.When the shear box started to move, shear strength was generated, which gradually increased and started to decrease after reaching a maximum value, and the simulation recorded the shear strength under different loadings.Figure 16 shows a comparison of the simulation and test results for the maximum shear strength under different vertical loadings.As the vertical loading increased, the maximum shear strength of the soil particles gradually increased, and the simulation and test results showed the same trend.When the vertical loading was 100 kPa, the average value of the simulation result was slightly higher than the test result, with a difference of 1.22 kPa.When the vertical loading was 200 kPa, the average value of the simulation result was also slightly higher than the test result, with a difference of 2.77 kPa.When the vertical loading was 300 kPa, the average value of the simulation result was slightly lower than the test result, with a difference of 1.26 kPa.When the vertical loading was 400 kPa, the average value of the simulation result was slightly higher than the test result, with a difference of 4 kPa.The analysis showed that the simulation and test values of the maximum shear strength of the soil particles were very close to each other for the same loading.The percentage difference between the simulation and the test results was greatest when the vertical loading was 200 kPa.The value was 8.8% < 10%, which is within the acceptable deviation range.
The results of the direct shear test were compared with the relevant reference [23,24].It could be seen that the maximum shear strength was different under different vertical loadings due to the different soil materials tested.However, the trends of the test results in this paper and the corresponding reference results were the same.Figure 16 shows a comparison of the simulation and test results for the maximum shear strength under different vertical loadings.As the vertical loading increased, the maximum shear strength of the soil particles gradually increased, and the simulation and test results showed the same trend.When the vertical loading was 100 kPa, the average value of the simulation result was slightly higher than the test result, with a difference of 1.22 kPa.When the vertical loading was 200 kPa, the average value of the simulation result was also slightly higher than the test result, with a difference of 2.77 kPa.When the vertical loading was 300 kPa, the average value of the simulation result was slightly lower than the test result, with a difference of 1.26 kPa.When the vertical loading was 400 kPa, the average value of the simulation result was slightly higher than the test result, with a difference of 4 kPa.The analysis showed that the simulation and test values of the maximum shear strength of the soil particles were very close to each other for the same loading.The percentage difference between the simulation and the test results was greatest when the vertical loading was 200 kPa.The value was 8.8% < 10%, which is within the acceptable deviation range.The above analysis showed that the results of the calibration parameters for the soil particles in this paper were correct.The parameters were taken as e = 0.3, λ = 0.35, n = 1, k = 0.67, μ1 = 0.9, μ2 = 0.7, and γ = 1 J/m 2 .The results of the direct shear test were compared with the relevant reference [23,24].It could be seen that the maximum shear strength was different under different vertical loadings due to the different soil materials tested.However, the trends of the test results in this paper and the corresponding reference results were the same.
The above analysis showed that the results of the calibration parameters for the soil particles in this paper were correct.The parameters were taken as e = 0.3, λ = 0.35, n = 1, k = 0.67, µ 1 = 0.9, µ 2 = 0.7, and γ = 1 J/m 2 .

Conclusions
In this paper, a systematic analysis of the test soil was carried out, and a geometric model of the soil was established.The EEPA model was identified as the contact model, and the parameters were calibrated.The accuracy of the calibrated parameters and the geometric model of the soil particles was verified.The main conclusions are presented below: (1) The geometry of the test soil was analyzed using an image particle analyzer, which approximated the soil particles as sphere-like and triangle-like and established a geometric model of the soil particles.
(2) Soil with a moisture content of 25% ± 1% was configured according to the test requirements; the density of soil particles was measured using the pycnometer method to be 1.844 g/cm 3 .
(3) The physical properties-compressibility and stickiness-of the soil used for the test were demonstrated by texture tests, and the EEPA model was chosen as the contact model for the soil particle simulation.
(4) A direct shear test was used to calibrate the corresponding parameters using shear strength as the response index.Firstly, the sensitivity of the parameters was analyzed through PB tests, and it was determined that the coefficient of static friction had a significant effect on the response index and that the surface energy and the coefficient of rolling friction had an approximately significant effect on the response index.Next, the BBD test was developed to calibrate the three parameters mentioned above by the response surface method, and the optimized parameter values were obtained.The applicability and accuracy of the parameters were further determined by comparing the simulation and test results of the direct shear test under different loadings, and the applicability of the soil particle model was also verified.

15 Figure 6 .
Figure 6.Soil moisture content and particle density test device.

Figure 6 .
Figure 6.Soil moisture content and particle density test device.

Figure 6 .
Figure 6.Soil moisture content and particle density test device.

Figure 8 .
Figure 8.The relationship between shear strength and displacement under different vertical loadings.

Figure 9 .
Figure 9.The relationship between the shear strength of the soil and the vertical loading.

Figure 8 .
Figure 8.The relationship between shear strength and displacement under different vertical loadings.

Figure 8 .
Figure 8.The relationship between shear strength and displacement under different vertical loadings.

Figure 9 .
Figure 9.The relationship between the shear strength of the soil and the vertical loading.

Figure 9 .
Figure 9.The relationship between the shear strength of the soil and the vertical loading.

Agriculture 2022 , 15 Figure 10 .
Figure 10.Diagram of direct shear box parts.The loading in the vertical direction was constant during the test, and the requirement was achieved through the secondary development of the API.The program flow chart is shown in Figure 11.In the diagram, F

Figure 10 .
Figure 10.Diagram of direct shear box parts.

Figure 10 .
Figure 10.Diagram of direct shear box parts.The loading in the vertical direction was constant during the test, and the requirement was achieved through the secondary development of the API.The program flow chart is shown in Figure 11.In the diagram, F

Figure 11 .
Figure 11.Flow chart of shear box working procedure.

Figure 11 .
Figure 11.Flow chart of shear box working procedure.

Figure 13 .
Figure 13.The relationship between force and displacement of the probe.

Figure 14 .
Figure 14.Pareto diagram of standardization effect of PBD test.

Figure 15 .
Figure 15.The relationship between the shear strength and the simulation time of the soil particles under different vertical loadings.

Figure 15 .
Figure 15.The relationship between the shear strength and the simulation time of the soil particles under different vertical loadings.

Agriculture 2022 , 15 Figure 16 .
Figure 16.Comparison of simulation and test results for shear strength under different vertical loadings.

Figure 16 .
Figure 16.Comparison of simulation and test results for shear strength under different vertical loadings.

Table 1 .
Mass distribution of soil particles in different particle size ranges.

Table 1 .
Mass distribution of soil particles in different particle size ranges.

Table 3 .
Analysis of variance of PBD test.
Figure 14.Pareto diagram of standardization effect of PBD test.

Table 5 .
Analysis of variance for BBD test.

Table 6 .
Analysis of model for BBD test.