Calibration and Testing of Parameters for the Discrete Element Simulation of Soil Particles in Paddy Fields

Abstract

The parameters of the discrete element simulation model for rice field soils serve as valuable data references for investigating the dynamic characteristics of the walking wheel of high-speed precision seeding machinery in paddy fields. The research specifically targets clay loam soil from a paddy field in South China. Calibration of essential soil parameters was achieved using EDEM_2022 software (and subsequent versions) discrete element simulation software, employing the Edinburgh Elasto-Plastic Adhesion (EEPA) nonlinear elastic-plastic contact model. The tillage layer and plough sub-base layer underwent calibration through slump and uniaxial compression tests, respectively. Influential contact parameters affecting slump and axial pressure were identified through a Plackett–Burman test. The optimal contact parameter combinations for the discrete element model of the tillage layer and plough sub-base layer were determined via a quadratic rotational orthogonal test. The accuracy of the discrete element simulation model’s parameters for paddy field soils was further validated through a comparative analysis of the simulation test’s cone penetration and the field soil trench test. Results indicate that the Coefficient of Restitution, surface energy, Contact Plasticity Ratio, and Tensile Exp significantly influence slump (p < 0.05). Additionally, the Coefficient of Restitution, Contact Plasticity Ratio, coefficient of rolling friction, and Tangential Stiff Multiplier significantly impact axial pressure (p < 0.05). Optimal contact parameters for the plough layer were achieved with a particle recovery coefficient of 0.49, a surface energy of 18.52 J/m², a plastic deformation ratio of 0.45, and a tensile strength of 3.74. For the plough subsoil layer, optimal contact parameters were a particle recovery coefficient of 0.47, a coefficient of interparticle kinetic friction of 0.32, a plastic deformation ratio of 0.49, and a tangential stiffness factor of 0.31. Results from the cone penetration test reveal no significant disparity in compactness between the actual experiment and the simulation test. The calibrated discrete element model’s contact parameters have been verified as accurate and reliable. The findings of this study offer valuable data references for understanding the dynamic characteristics of the walking wheel of the entire machinery in high-speed precision seeding in paddy fields.

1. Introduction

Rice cultivation constitutes approximately 25% of the total cultivated land dedicated to food crops in this nation [1]. The mechanization of paddy fields plays a pivotal role in advancing the modernization of our agriculture [2]. Nevertheless, current operations of paddy field machinery encounter challenges due to the resistance generated by mud and the interaction between the driving wheels and the paddy field plough subsoil. These factors result in soil pressure and shearing forces, potentially causing soil deformation and failure. The muddy terrain presents a formidable obstacle for machinery, leading to issues such as wheel slippage, head warping, and wheel entrapment in the mud, among other failure phenomena.

In recent years, significant attention has been directed towards advancing agricultural modernization. Consequently, the discrete element method has been employed to optimize and analyze agricultural machinery and the interactions between particles in the law of motion [3]. This method proves particularly advantageous and has been instrumental in constructing a soil model for paddy fields, facilitating the design and optimization of technical aspects related to paddy field machinery operation [4,5,6]. The discrete element method is utilized to explore the interaction between paddy field operating machinery and soil. Ensuring the reliability of the model necessitates precise calibration of simulation parameters, with scholars recommending the implementation of both actual and simulation tests for accuracy [7,8,9,10,11].

Shi et al. [12] combined the benefits of the Hysteretic Spring Contact Model (HSCM) and the Linear Cohesion Model (LCM), utilizing the soil accumulation angle as the response value to construct a model of farmland in the Northwest Dry Zone. Zhang et al. [13], employing the Hertz-Mindlin (no slip) contact model, developed a discrete elemental model of sandy soil, calibrating parameters through the angle of stacking test. Li et al. [14] modeled clay-heavy black soil particles in the northeast region using the Hertz-Mindlin with JKR contact model, with the simulated accumulation angle of soil particles serving as the response surface. Xiang et al. [15] combined Hertz-Mindlin with JKR contact models, calibrated parameters through stacking tests, and obtained accurate contact parameters for the discrete element simulation model of southern clay loam soil. Zhou et al. [16], selecting the Hertz-Mindlin model with JKR contact, calibrated parameters through a soil slump test, and derived optimum contact parameters for discrete elements of soil with high moisture content using the steepest-climbing test and central composite test. Xie et al. [17] opted for the Edinburgh Elasto-Plastic Adhesion contact model, acquiring the best contact parameter combinations for discrete elements in viscous and plastic deformation-prone soil through the uniaxial confined compression test and the unconfined compressive strength test, subsequently constructing a discrete element model of viscoelastic soil.

In summary, we propose the Edinburgh Elasto-Plastic Adhesion (EEPA) contact model, incorporating soil elasticity and viscosity, building upon the work of Thakur et al. [18]. However, there is a notable lack of literature addressing the calibration of discrete elemental parameters for the integration of cultivated and ploughed subsoil layers in paddy fields in South China. To address this gap, our study employs a combination of physical tests and discrete element simulation tests to calibrate parameters for clay loam soil with high water content. The discrete element soil bed is established using the EEPA contact model, and the particle contact parameters of the model are refined by comparing actual and simulated results from the slump test and axial compression test. This calibration process leads to the development of a discrete element simulation model specifically tailored for paddy field soil. To validate parameter accuracy, cone penetration tests are conducted. Our objective is to create a preliminary discrete element model of clay loam soil in paddy fields, providing valuable data references for optimizing the structural design of lightweight power chassis in paddy fields.

2. Materials and Methods

2.1. Soil Basic Parameters

Soil samples were selectively extracted from the experimental plot situated at the teaching and research center of South China Agricultural University (SCAU). The analysis of the soil texture identified it as clay loam, characterized by elevated levels of porosity, water content, and plasticity. A total of nine distinct soil sample groups were randomly procured from the experimental field. The determination of moisture content, density, and liquid-plastic limit of the soil samples was conducted using the drying method, ring knife method, and liquid-plastic limit measurement, respectively. The outcomes of these analyses are depicted in Figure 1.

2.2. Contact Model Theory and Parameters to Be Calibrated

The clay present in paddy field planting exhibits robust adhesion and elasticity. The Edinburgh Elasto-Plastic Adhesion (EEPA) non-linear elastic-plastic contact model, employed in the EDEM software, aptly captures the characteristics of soil particles both before and after compaction. Post-compaction, the particles acquire a cohesive and deformable nature yet retain a discrete state when in an uncompacted condition [19].

The normal spring force of the EEPA contact model, seamlessly integrated into the EDEM_2022 software (and subsequent versions), is illustrated in Figure 2 and mathematically expressed by the following equation [18].

The comprehensive standard normal force $f_n$ results from the combination of the hysteresis elastic force $f_{hys}$ and the normal damping force $f_{nd}$. The particle overlap relationship governing the normal contact force is expressed as follows [18]:

$$\begin{array}{c}{f}_n=\left(f_{h{y}_s}+f_{nd}\right)\end{array}$$

(1)

$$\begin{array}{c}{f}_{hys}=\left\{\begin{array}{c}{f}_0+k_1{\delta }^n if k_2\left({\delta }^n{-\delta }_p^n\right)\ge k_1{\delta }^n\\ f_0+k_2\left({\delta }^n-{\delta }_p^n\right) if k_1{\delta }^n> k_2\left({\delta }^n-{\delta }_p^n\right)>-k_{adh}{\delta }^n\\ f_0-k_{adh}{\delta }^n if-k_{adh}{\delta }^n\ge k_2\left({\delta }^n-{\delta }_p^n\right)\end{array}\right.\end{array}$$

(2)

where $u$ represents the unit normal vector from the contact point to the center of mass, $f_0$ denotes the initial bond strength of the particles, $k_1$ refers to the loading stiffness coefficient, $k_2$ is the unloading stiffness coefficient, ${\delta }_p$ represents the contact overlap of the particles, and $k_{adh}$ is the coefficient of viscous strength.

$$\begin{array}{c}{f}_{nd}={\beta }_n{v}_n\end{array}$$

(3)

$$\begin{array}{c}{\beta }_n=\sqrt{\frac{4m^{\ast }{k}_1}{1+{\left(\frac{\pi }{\mathit{ln}e}\right)}^2}}\end{array}$$

(4)

$$\begin{array}{c}{m}^{\ast }=m_i·m_j/\left(m_i+m_j\right)\end{array}$$

(5)

where the normal velocity is represented by $v_n$, the damping coefficient by ${\beta }_n$, the contact particle by $m^{\ast }$, and the coefficient of recovery by $e$.

Tangential damping is calculated by multiplying the interparticle tangential damping coefficient ${\beta }_t$ with the relative tangential velocity $v_t$. The t description of the tangential damping coefficient can be found in reference [18].

$$\begin{array}{c}{\beta }_t=\sqrt{\frac{4m^{\ast }{k}_t}{1+{\left(\frac{\pi }{\mathit{ln}e}\right)}^2}}\end{array}$$

(6)

The computation of the maximum tangential friction adheres to the Coulomb friction criterion, with adjustments made to the tangential normal force based on the bonding force specified in this equation:

$$\begin{array}{c}{f}_{ct}\le u\left(\left|f_{hys}+k_{adh}{\delta }^n-f_0\right|\right)\end{array}$$

(7)

where $f_{ct}$ represents the maximum tangential friction and $u$ denotes the friction coefficient.

In summary, when applying the EEPA contact model, the requisite input parameters include the Coefficient of Recovery (e), coefficient of static friction (μ_s), coefficient of rolling friction ($\mu$_r), Constant pull-off force ($f_0$), surface energy ($\mathsf{\Delta }\gamma$), Contact Plasticity Ratio (λ_ρ), Tensile Exp (R_m), Slope Exp (n), and Tangential Stiff Multiplier (K_tm). It is crucial to utilize this set of parameters for the accurate implementation of the EEPA contact model. To streamline calibration, the value of $f_0$ is set to 0. The loading branch index can be either 1 or 1.5, with the latter being optional for a more precise representation of the nonlinear stress-strain properties of the soil. Therefore, the selection of n is recommended to be 1.5 [19].

2.3. Characteristics of Soil Layers in Paddy Fields

2.3.1. Tillage Calibration Test

The angle of repose test is a common method for calibrating model parameters for discrete elements [20,21,22,23]; However, its applicability is limited in soils with high water content. Consequently, the slump test was chosen for its ability to provide a reliable representation of the soil’s cohesion, friction properties, and mobility. The testing apparatus, as illustrated in Figure 3, utilizes a slump barrel tester with a height of 300 mm, an upper diameter of 100 mm, and a lower diameter of 200 mm. Tillage layer soil samples are placed into the standard cone slump barrel according to the prescribed method and thoroughly scraped. The barrel is then lifted vertically, and the soil samples naturally slump due to gravity. The downward slumping is measured, with this dimension (mm) serving as the slump, indicating fluidity, where a larger slump signifies better fluidity.

The simulation parameters for the slump test are detailed in

Table 1. Simulation test parameters.

Materials	Parameter	Unit	Value	Source
Tillage soil-soil	Poisson’ s Ratio (ν)		0.3	[9]
	Solids Density ($\rho$)	g/cm3	1.7	measurement
	Shear Modulus (G)	Pa	5.0 × 105	[9]
Plow subsoil—soil	ν		0.38	[24]
	$\rho$	g/cm3	2.2	measurement
	G	Pa	1.3 × 106	measurement
EquipMlaterial	ν		0.3	[25]
	$\rho$	g/cm3	1.2	[25]
	G	Pa	3.16 × 109	[25]
Tillage soil-soil	Coefficient of Restitution e			To be calibrated
	Coefficient of static friction (μs)			To be calibrated
	Coefficient of rolling friction (μr)			To be calibrated
Plow subsoil—soil	e			To be calibrated
	μs			To be calibrated
	μr			To be calibrated
Soil-EquipMlaterial	e		0.3	[15]
	μs		0.6	[15]
	μr		0.1	[15]
Contact model	Surface energy (∆γ)			To be calibrated
	Contact Plasticity Ratio ($\lambda \rho$)			To be calibrated
	Tensile Exp (Rm)			To be calibrated
	Tangential Stiff Multiplier (Ktm)			To be calibrated

, with the EEPA contact model chosen for particle-particle contact. The particle-contact component is selected as the Hertz-Mindlin (no slip) contact model to simplify the simulation, with the contact component and particle viscosity factor excluded. The simulation generates particles to fill a cylinder that is 300 mm in height and has an upper diameter of 100 mm and a lower diameter of 200 mm. The cylinder is raised at a rate of 30 mm/s. The simulation concludes once the particles cease movement, and the slump is determined by measuring the cylinder’s height and calculating H = 300 − h.

2.3.2. Calibration Test for Plough Base Layer

As the machine traverses the paddy field, it compacts the soil particles, inducing plastic deformation and generating reaction forces toward the wheels. When subjected to the uniaxial confined compression test, the soil undergoes elastic-plastic deformations, resulting in an axial stress-strain curve that reflects the load-bearing characteristics of the subsoil under the plough.

In the experiment, a soil specimen was placed into a 45 mm diameter cylinder and tested using a WD-E precision micro-control electronic universal testing machine. The pressure disc was connected to the pressure meter of the universal tester and subjected to a constant downward load of 10 mm/s. The load pressure was recorded when the downward compression depth reached 8 mm. Subsequently, the original sample was extracted, and the axial stress-strain chart was generated using the universal testing machine (Figure 4a).

For the axial compression simulation test, the simulation parameters are outlined in Table 1. A cylindrical model with a height of 250 mm and a bottom diameter of 45 mm (adjusted from the actual testing cylinder height of 100 mm) is simulated, as depicted in Figure 4b. The loading rate is maintained at 10 mm/s, consistent with the actual rate. The platen halts downward motion when the cylinder particles are compressed to 102 mm, and the axial stress-strain curves are obtained through the post-processing function of EDEM.

2.4. Significance of DEM Parameters on Particle Response

Considering the characteristics of the EEPA contact model, the associated parameters underwent screening via a Plackett–Burman test. Viscous and elastic-plastic indicators were employed as response values during the simulation test. The levels of each of the seven factors were established, as outlined in

Table 2. Parameters of the Plackett–Burman test.

Parameter	Factors	Levels
		−1	1
X1	Coefficient of Restitution	0.2	0.6
X2	coefficient of static friction	0.3	0.9
X3	coefficient of rolling friction	0.2	0.6
X4	surface energy	10	30
X5	Contact Plasticity Ratio	0.2	0.6
X6	Tensile Exp	1.5	4
X7	Tangential Stiff Multiplier	0.3	0.9
X8–X11	Virtual parameter	——	——

. Additionally, four virtual parameters were included for blank control error analysis, with each parameter assuming both high and low levels within its range of values.

2.5. Center Composite Design Experiment

The Central Composite Design (CCD) employs a 2-level factorial design with extra centroids and axial points to fit a quadratic model. The conventional CCD typically comprises five levels for each factor, but this can be adjusted by selecting an axial distance of 1.0, resulting in a centered composite design on the face with only three levels per factor. Reproducing the centroid enhances the accuracy of predictive capabilities in the vicinity of the factor space’s center. CCD experiments were conducted with the simulation parameters to anticipate the correlation between the particles’ physical features and the key factors influencing the Discrete Element Method (DEM) parameters. The simulation setup remained consistent with the description in Section 2.3, and the coding tables for the compression test factor levels are presented in

Table 3. Factor level coding table for the quadratic orthogonal rotary combination test of the tillage layer.

Levels	X1	X4	X5	X6
γ	0.68	34.14	0.68	4.5
1	0.6	30	0.6	4
0	0.41	20	0.41	2.75
−1	0.22	10	0.22	1.5
−γ	0.14	5.86	0.14	0.98

, while the coding table for the axial compression test factor levels is shown in

Table 4. Factor level coding table for the quadratic orthogonal rotary combination test of the plow pan.

Levels	X1	X3	X5	X7
γ	0.68	0.68	0.68	0.95
1	0.6	0.6	0.6	0.9
0	0.41	0.41	0.41	0.6
−1	0.22	0.22	0.22	0.3
−γ	0.14	0.14	0.14	0.26

.

3. Results

3.1. Measurements of Soil Physical Parameters

The settlement test was replicated three times in Section 2.3.1 with a settlement (H) of 140 mm, as illustrated in Figure 3a. Additionally, in Section 2.3.2, the axial pressure was recorded at 180.24 N when the soil sample was compressed to 8 mm, as depicted in Figure 4b. The outcomes of these two sets of tests serve as the target response values calibrated to the simulation parameters.

3.2. Significance of the Effect of DEM Parameters on Particles

The results of the Plackett–Burman test were analyzed using Design-Expert_10 software to generate Pareto plots and elucidate the impact of Discrete Element Method (DEM) parameters on the response values. A DEM parameter is deemed to have a highly significant effect on the response value when its contribution surpasses the Bonferroni limit. Conversely, it has almost no effect when its contribution falls below the t-value limit.

Figure 5 illustrates that the Coefficient of Restitution, Contact Plasticity Ratio, coefficient of rolling friction, and Tangential Stiff Multiplier exhibit significance for axial compression. Additionally, the Contact Plasticity Ratio is highly significant for deflection, while the Coefficient of Restitution, surface energy, Contact Plasticity Ratio, and Tensile Exp are significant for slump. The coefficient of static friction has a minimal effect on both axial pressure and deflection. Consequently, the Coefficient of Restitution, Contact Plasticity Ratio, coefficient of rolling friction, and Tangential Stiff Multiplier were utilized as independent parameter variables to predict axial pressure variation. Moreover, the Coefficient of Restitution, surface energy, Contact Plasticity Ratio, and Tensile Exp were used as single parameter variables to predict changes in slump.

3.3. Calibration Parameters

3.3.1. Tillage Parameters

Quadratic multiple regression analysis, conducted with Design Expert software, was employed to generate response surfaces and contour plots (Figure 6) depicting the interactions among parameters influencing the sink. In Figure 6a, the impact of the X₁–X₄ interaction on the sink is illustrated. The response surface curve for X₁ is steeper than that of X₄, and the contour density along X₁ is slightly higher than the density along X₄. This observation signifies that X₁ has a more substantial effect on the sink.

For the experimental results presented in

Table 5. Design and results of the quadratic orthogonal rotational combination test for the slump test.

No.	X1	X4	X5	X6	Slump (mm)
1	0.22	10	0.22	1.5	195
2	0.6	10	0.22	1.5	185
3	0.22	10	0.6	1.5	190
4	0.6	10	0.6	1.5	155
5	0.22	30	0.22	1.5	200
6	0.6	30	0.22	1.5	180
7	0.22	30	0.6	1.5	230
8	0.6	30	0.6	1.5	40
9	0.22	10	0.22	4	200
10	0.6	10	0.22	4	195
11	0.22	10	0.6	4	160
12	0.6	10	0.6	4	145
13	0.22	30	0.22	4	205
14	0.6	30	0.22	4	195
15	0.22	30	0.6	4	215
16	0.6	30	0.6	4	50
17	0.14	20	0.41	4	170
18	0.68	20	0.41	2.75	145
19	0.41	20	0.14	2.75	210
20	0.41	20	0.68	2.75	160
21	0.41	5.86	0.41	2.75	190
22	0.41	34.14	0.41	2.75	147
23	0.41	20	0.41	1	153
24	0.41	20	0.41	4.5	155
25	0.41	20	0.41	2.75	150
26	0.41	20	0.41	2.75	152
27	0.41	20	0.41	2.75	148

, the regression equation was retained using Design Expert software under the condition of non-significance (p > 0.05) in the Analysis of Variance (ANOVA):

$$\begin{array}{l}S=134.61159+220.47096X_1+2.60565X_4-17.56948X_5+10.48493X_6\\ -10.52633X_1{X}_4-623.2684X_1{X}_5+21.09022X_1{X}_6-3.94737X_4{X}_5+0.2X_4{X}_6\\ -21.05263X_5{X}_6+67.54763X_1^2\end{array}$$

(8)

where S is the slump.

By solving Equation (8) with the measured S as the objective, 73 sets of solutions were obtained. These solutions are essentially indistinguishable from each other. Utilizing the desirability values as a reference, a set of parameters X₁ = 0.49, X₄ = 18.52 J/m², X₅ = 0.45, and X₆ = 3.74 is selected, as it has the desirability values closest to 1, indicating the highest degree of reliability.

3.3.2. Parameters of the Plough Substrate

Quadratic multiple regression analysis, facilitated by Design Expert software, was employed to create response surfaces and contour plots illustrating the interactions among parameters influencing slump (Figure 7). For the experimental results detailed in

Table 6. Design and results of quadratic orthogonal rotational combination tests for axial force tests.

No.	X1	X3	X5	X7	Axial Force (N)
1	0.22	0.22	0.22	0.3	246.5
2	0.22	0.22	0.22	0.3	249.7
3	0.22	0.22	0.6	0.3	187.3
4	0.6	0.22	0.6	0.3	179.5
5	0.22	0.6	0.22	0.3	234.1
6	0.6	0.6	0.22	0.3	216.6
7	0.22	0.6	0.6	0.3	174.12
8	0.6	0.6	0.6	0.3	143.8
9	0.22	0.22	0.22	0.8	235.6
10	0.6	0.22	0.22	0.8	223.95
11	0.22	0.22	0.6	0.8	146.3
12	0.6	0.22	0.6	0.8	135.7
13	0.22	0.6	0.22	0.8	240.4
14	0.6	0.6	0.22	0.8	212
15	0.22	0.6	0.6	0.8	169.5
16	0.6	0.6	0.6	0.8	126.36
17	0.14	0.41	0.41	0.55	223.13
18	0.68	0.41	0.41	0.55	190.5
19	0.41	0.41	0.14	0.55	225.2
20	0.41	0.41	0.68	0.55	113.7
21	0.41	0.14	0.41	0.55	168.5
22	0.41	0.68	0.41	0.55	154.68
23	0.41	0.41	0.41	0.2	176.1
24	0.41	0.41	0.41	0.9	153.9
25	0.41	0.41	0.41	0.55	155.2
26	0.41	0.41	0.41	0.55	153.9
27	0.41	0.41	0.41	0.55	155.3

, the regression equation was retained using Design Expert software under the condition of non-significance (p > 0.05) in the Analysis of Variance (ANOVA):

$$\begin{array}{l}F=440.9886-498.1707X_1-293.14681X_3-120.85088X_5-116.49822X_7\\ -70.34421X_1{X}_3-165.56581X_1{X}_5-58.54055X_1{X}_7+36.90254X_3{X}_5\\ -89.{76986X}_3{X}_7-137.86172X_5{X}_7+704.06121X_1^2+704.06121X_3^2\\ +83.69084X_5^2+77.56852X_7^2\end{array}$$

(9)

where F is the axial pressure.

Several sets of solutions are obtained by solving Equation (9) with the measured F objective. These solutions are almost indistinguishable from each other, and using the desirability values as a reference, a set of parameters X₁ = 0.47, X₃ = 0.48, X₅ = 0.32, and X₇ = 0.31 is chosen, as it possesses desirability values closest to 1, signifying the highest reliability.

3.4. Simulation Verification of an Optimal Parameter Set

Discrete element simulation modeling was conducted using EDEM2021 software, leveraging the outcomes of cohesion and elastic-plastic parameter calibration tests. To assess the precision of the calibration results, slump and axial compression tests were executed. The slump measured 143 mm, with a minimal error of only 2.14%, while the axial compression registered 175.613 N, with an error of merely 2.57%. These results suggest that a combination of significance analysis and the response surface method can effectively optimize the physical parameters of soil particle simulation.

3.5. Simulation Results and Discussion

3.5.1. DEM Modelling

The soil layer in the paddy field is primarily divided into three layers, necessitating the establishment of a three-layer DEM model. The tillage layer and the plough subsoil layer directly interact with agricultural equipment and carriers. The tillage layer has a thickness of 180–220 mm, the plough subsoil layer ranges from 60–100 mm, and the combined thickness of the two layers is 240–320 mm. For ease of observing mechanical-particle interactions, minimizing soil particle extrusion by the soil trough shell, and optimizing computational efficiency, the soil trough dimensions were set to length: 3200 mm and width: 1800 mm, as depicted in Figure 8. This configuration is designed to establish a DEM model suitable for clay loam soils in South China.

3.5.2. Testing for Model Validation

To further validate that the discrete element particle model, constructed after parameter optimization, accurately reflects the physico-mechanical parameters of actual field soils, the cone penetration test is employed. Cone penetration is a highly intuitive indicator of the overall mechanical properties of the soil. This test, in conjunction with soil trench field tests and discrete element simulation tests [26], enables a comprehensive assessment. Therefore, the relative error value between the measured value and the simulation result in the cone penetration test is utilized to judge the accuracy of the discrete element particle model (Figure 9).

Compaction tests were conducted to validate the mechanical accuracy of the constructed discrete element model. A soil stiffness meter was utilized to measure the stiffness at various soil depths (0–40 cm), with data collected from five points in the test field and averaged. In EDEM, a 1:1 soil compactometer probe model was imported, and the probe was positioned on the surface of the soil particles with a descent rate of 30 mm/s². At the conclusion of the simulation test, the total contact force between the probe model and the soil particles in the vertical direction was exported to the post-processing module. The comparison results revealed no significant difference in tightness, affirming the suitability of this DEM model for the test field.

4. Discussion

(1) The static and rolling friction coefficients between soil particles play a pivotal role in shaping their behavior. Wang et al. [27] established a method for acquiring the coefficient of static friction between particles and plates using a platform. However, determining the static and rolling friction coefficients in particle-particle interactions through experimental means poses challenges. To enhance simulation accuracy and alleviate computational demands, a common strategy is to select a subset of crucial parameters. Su et al. [28] found that the combination of factors is a significant parameter for designing experiments using Plackett–Burman.

Initially, we employed the Plackett–Burman experiment to identify the noteworthy influencing factors. Subsequently, we determined specific parameter values using the Central Composite Design (CCD) method. The study underscores that the rolling friction coefficient between soil particles significantly impacts axial pressure, while none of the static friction coefficients exert a notable effect. The parameters governing soil particle interactions were discerned through a blend of experimental measurements and simulation optimization.

(2) The angle of repose test is a widely employed method for calibrating discrete elemental model parameters [20,21,22,23]. Nevertheless, its applicability is constrained when dealing with soils characterized by high water content, primarily due to inherent limitations. Zhou’s [16] findings revealed that the slump test is an accurate calibration method for soils with high water content. Consequently, we selected this test method for calibrating soils with high water content in the till layer. However, in our approach, we chose the EEPA contact model over the JKR contact model. The discrete element slump simulation, calibrated through the Plackett–Burman test and the Central Composite Design (CCD), exhibited a relative error of 2.14% compared to the actual test results. This confirms the suitability of the obtained parameters for constructing the discrete element model. To ensure accurate representation of soil particles, validation of the physical and mechanical parameters of the model was conducted through cone penetration testing, performed subsequent to the establishment of the DEM simulation model.

5. Conclusions

In this study, the parameters of the soil layer EEPA model were calibrated using the Plackett–Burman test and the quadratic rotated orthogonal test. The results of the calibration and validation tests were analyzed, leading to the following conclusions:

(1)

The results of the Plackett–Burman test showed that the Coefficient of Restitution, Contact Plasticity Ratio, coefficient of rolling friction, and Tangential Stiff Multiplier in the contact model were significant for axial pressure. Contact Plasticity Ratio was highly significant for deflection. The Coefficient of Restitution, surface energy, Contact Plasticity Ratio, and Tensile Exp were significant for deflection. The coefficient of static friction had almost no effect on axial pressure and slump.

(2)

The regression model was solved and fitted to the measured values to give the Coefficient of Restitution, Contact Plasticity Ratio, Tensile Exp, and surface energy of 0.49, 0.45, 3.74, and 18.52 J/m², respectively; The Coefficient of Restitution, Contact Plasticity Ratio, coefficient of rolling friction, and Tangential Stiff Multiplier were 0.47, 0.49, 0.32, and 0.31, respectively. The slump obtained through constructing the discrete element soil model with the optimal parameter set exhibits a relative error of 2.14% compared to the measured value, while the axial pressure demonstrates a relative error of 2.57% from the measured value.

(3)

The optimized discrete element model underwent cone penetration tests and validation through field trials. The results of the validation test revealed no significant difference in soil particle compactness, indicating that the model can accurately simulate soil mechanical properties. This model serves as a valuable reference for further research on the dynamic characteristics of the entire traveling wheel of high-speed precision seeding machinery in paddy fields.

This study makes a significant contribution to the ongoing research on constructing a soil model for paddy fields and analyzing the kinetic characteristics of the entire walking wheel of high-speed precision seeding machinery in paddy fields. Its specific emphasis lies in investigating the mechanical characteristics and physical traits of soil particles during the operational process of agricultural machinery. However, the primary focus of this paper is on presenting a method for constructing a discrete element particle contact model for soil and calibrating its parameters. It’s noteworthy that the model parameters exhibit complexity and variability across different soil properties and water content conditions in various regions. Therefore, reconstruction and calibration become essential, considering factors such as soil type, particle size composition, and water content.