The Calibration of Soil Simulation Parameters for Wheat Grown after Rice in the Yangtze River Basin of China

Abstract

In recent years, there have been many studies on the calibration of soil simulation parameters; however, there are few soil parameters for wheat grown after rice that have been calibrated in the lower reaches of the Yangtze River, and the data from such calibrations remain inaccurate. Therefore, using the soil of Jiangsu as our research object, a soil parameter calibration was carried out based on the EEPA (Edinburgh elastoplastic adhesion) model and using the discrete element software EDEM (2020, DEM-Solutions, Edinburgh, UK). The depression depth measured via a uniaxial compression test and the maximum crushing force measured via an unconfined compression test were taken as indexes. The Plackett–Burman test was used to screen seven influencing factors (recovery coefficient, static friction coefficient, rolling friction coefficient, surface energy, contact plastic deformation ratio, tensile exp, and tangential stiff multiplier). The recovery coefficient and static friction coefficient were analyzed using a central composite test with depression depth as an index. The surface energy, plastic deformation ratio, and tangential stiffness factor were analyzed via a Box–Behnken test, with the maximum crushing force as the index. Taking the measured depression depth of 3.36mm and the maximum destructive power of 210 N as the target, the following final data were obtained—recovery coefficient: 0.322, static friction coefficient: 0.676, rolling friction coefficient: 0.5, surface energy: 17.158, contact plastic deformation ratio: 0.358, tensile exp: 2, and tangential stiff multiplier: 0.8. Finally, the simulation value and the actual value of the parameter group were verified and compared. It was found that the mismatch coefficient R² of the actual value and the simulation value is 93.509%. The mismatch coefficient R² between the actual and simulated values of unconfined compressive test is 94.2%. This shows that the curves obtained from the real test and simulation test have a high similarity. This study can provide technical support for the simulation and optimization of growing wheat after rice seeding equipment in the lower reaches of Yangtze River in China.

1. Introduction

Wheat is an important food crop in China as it contributes to food security [1]. As a major agricultural province, Jiangsu’s main cropping system is a rice–wheat rotation, with rice grown after wheat accounting for over 70% of the total area initially used for growing wheat [2]. Jiangsu is also a large contributor to the production of agricultural machinery. In 2023, the comprehensive mechanization rate of Jiangsu’s crop tillage and harvesting reached more than 87 percent, and wheat yields reached about 13.5 million tons. Although mechanical utilization rates and wheat yield are high, there are still some challenges faced in the sowing stage, such as hanging and drying seeds. Therefore, sowing and seeding equipment must be optimized urgently to increase wheat yield. However, the optimization of machinery requires field testing, which is expensive, time-consuming, and has many limitations depending on the weather and soil conditions. Many experts and scholars use simulation methods to address these shortcomings [3]. However, most of the soil simulated by domestic experts is for Hebei, Hunan, and other regions of China. There are few studies on the calibration of soil simulation parameters for wheat in the lower reaches of the Yangtze River, and the parameters of different soils are different. The accurate calibration of simulation parameters is the key to a successful simulation. Therefore, we conducted a calibration of simulation parameters for rice–stubble wheat soil in Jiangsu, thus enriching the soil parameters of rice–stubble wheat in the Yangtze River Basin of China. The data that were found support the optimization of the rice–stubble wheat seeding mechanism.

In recent years, scholars at home and abroad have used the discrete element simulation method for the detailed study of sowing operation technology, determining the interaction between soil and machinery [4,5,6,7,8] to better optimize machinery [9,10]. Wang Xianliang et al. [11] argue that a discrete element soil model can be established based on the parameters measured in basic experiments. Combined with the stacking angle and shear tests, a model simulation can be used to analyze the sensitivity of the model’s parameters. With sensitivity parameters as variables and real test measurements as target values, a proxy model was constructed, and parameters were optimized using a Gauss–Newton iteration. The results show that the calibrated optimized soil model can effectively replace real soil in simulations. Through a uniaxial closed compression test and an unconfined compressive strength test, Xie Fangping et al. [12], taking unconfined compressive strength as the viscosity index and axial strain as the plastic index, conducted two central composite tests, established a quadratic regression model, optimized the solution, and obtained the optimal parameter combination. Zhan Zhao et al. [13] applied the EEPA (Edinburgh elastoplastic adhesion) model to analyze the relationship between compressive strength under different advancing speeds and compaction depths and the soil compaction degree, establishing a relationship model between compressive strength under different advancing speeds and the compaction depth and the soil compaction degree. Ultimately, according to the set compression strength, the nursery bed compaction degree can be automatically adjusted in real time. Wanli Pengcheng et al. [14] used the EEPA model to establish the dimensions of four kinds of soil particles—single-sphere, double-sphere, three-sphere, and four-sphere—as well as six kinds of virtual soil tanks by using the following filling methods: a single-sphere, single-sphere variable particle size, and a double-sphere, three-sphere, four-sphere and multi-sphere mixture. By measuring the static stacking angle, axial pressure, tillage resistance, and tilling stacking angle of the soil in the simulated soil tank and comparing them with those of the indoor test results, it was concluded that the virtual soil tank should be constructed using the four-sphere or multi-sphere mixed filling method in the simulation of the tilling accumulation angle. Kim et al. [3] used the discrete element method to model agricultural soil, measured soil properties with tillage depth (0–200 mm) as the target, and used the virtual leaf shear test and cone penetration test values, instead of the resting angle test, to calibrate the DEM soil model. Finally, field tests were carried out to verify the accuracy of the tractive force prediction simulation.

This paper takes soil for growing wheat after rice in Jiangsu as the research object, uses the EEPA model, and addresses the subsidence depth and maximum crushing force. The Plackett–Burman, central composite, and Box–Behnken tests were used to carry out parameter calibration for seven influencing factors (recovery coefficient, static friction coefficient, rolling friction coefficient, surface energy, contact plastic deformation ratio, tensile exp, and tangential stiff multiplier). The comparison and verification with the actual values provide the basis for the parameters, a feasible calibration method for the discrete element simulation research based on the EEPA model, and technical means for the optimization of machine tools.

2. Materials and Methods

It is difficult to measure soil particle parameters directly in the process of establishing a soil model via the discrete element method (DEM). If the discrete element soil model is established based on indirectly measured soil parameter values, the error in the simulation results will be large. Therefore, experts and scholars often use calibration methods to measure soil parameters. In recent years, the number of DEM-related studies has rapidly increased, and new methods of parameter calibration are continuously emerging [15]. Contact models commonly used in soil calibration include the Hertz–Midlin (no slip), JKR (Johnson–Kendall–Roberts), bonding, hysteretic spring + linear cohesion, and EEPA models. The Hertz–Midlin (no slip) model is mainly used for particles without special forces between them and is an ideal state model [16]; the JKR model is mainly used for particles with adhesion between surfaces [17]; the bonding model is mainly used in the presence of fragmentation within particles, which is characterized by the failure of the bound bond to regenerate after fracture; the hysteretic spring + linear cohesion model is a combination of the elastoplastic and adhesive force, and the elastoplastic relationship is linear [6,18]; and the EEPA model combines the elastoplasticity of its own particles with the adhesion between particles. Its elastoplastic deformation is nonlinear, which is more consistent with the modeling of soil, especially compressible clayey soil [19]. The particle contact and normal force overlap curves of the EEPA model are shown in Figure 1 [20]. When two discrete soil particles are pressed together, they undergo elastic and plastic deformation, and the N-bond force increases with the increase in the size of the plastic contact area.

f_n is the normal bonding force, N; f₀ is the constant drag force, N; f_min is the maximum adhesion between particles, N; δ is the overlap length; δ_p is the plastic deformation of particles; δ_min is the amount of overlap of particles with maximum adhesive force; δ_max is the maximum length of the overlap area; k₁ is the original loading stiffness; K₂ is the reloading stiffness; k_adh is the attenuation stiffness; X is the adherent branch curve power index; and n is the nonlinear curve power index.

The total contact normal force f_n is the sum of the hysteresis spring force f_hys and the normal damping force f_nd, which can be written as

f_n = f_hys + f_nd

(1)

f_adh and f_nd can be derived from the following equation:

$${\mathrm{f}}_{\mathrm{h}\mathrm{y}\mathrm{s}}\left\{\begin{array}{c}{\mathrm{f}}_0+{\mathrm{k}}_1{\mathsf{\delta }}^{\mathrm{n}} {\mathrm{k}}_2\left({\mathsf{\delta }}^{\mathrm{n}}-{\mathsf{\delta }}_{\mathrm{p}}^{\mathrm{n}}\right)\\ {\mathrm{f}}_0+{\mathrm{k}}_2\left({\mathsf{\delta }}^{\mathrm{n}}-{\mathsf{\delta }}_{\mathrm{p}}^{\mathrm{n}}\right) {\mathrm{k}}_1{\mathsf{\delta }}^{\mathrm{n}}>{\mathrm{k}}_2\left({\mathsf{\delta }}^{\mathrm{n}}-{\mathsf{\delta }}_{\mathrm{p}}^{\mathrm{n}}\right) \\ {\mathrm{f}}_0-{\mathrm{k}}_{\mathrm{a}\mathrm{d}\mathrm{h}}{\mathsf{\delta }}^{\mathrm{x}} -{\mathrm{k}}_{\mathrm{a}\mathrm{d}\mathrm{h}}{\mathsf{\delta }}^{\mathrm{x}}\ge {\mathrm{k}}_2\left({\mathsf{\delta }}^{\mathrm{n}}-{\mathsf{\delta }}_{\mathrm{p}}^{\mathrm{n}}\right)\end{array}\right.$$

(2)

where δ is the contact overlap, f₀ is the constant tension, k₁ is the initial loading stiffness, k₂ is the unloading/reloading stiffness, k_adh is the bonding stiffness (k_adh = f_min/δ_minx), n is the power of k₁ and k₂, x is the power of k_adh, and v_nrel is the normal component of the relative velocity.

2.1. Materials

The soil test field is located in Baima Base, Nanjing Institute of Agricultural Mechanization, Lishui District, Nanjing, Jiangsu Province. This area is located in the lower reaches of the Yangtze River and has a subtropical monsoon climate. The average annual temperature is 15.8 degrees Celsius, with an average of 2008 hours of sunshine, and an annual precipitation of 1048.6 mm. Rice and wheat crops are rotated here all year round.

In this study, a ring knife with an inner diameter of 50 mm and a height of 15 mm was used to randomly sample the soil in this area [21]. The ring knife was greased before sampling to facilitate the extraction of excavated samples. The obtained soil samples (15 samples) were sealed, stored, and brought back to the laboratory for soil parameter testing.

2.2. Soil Intrinsic Parameter Measurement

The samples were cut into appropriate sizes and weighed to obtain the quality of the unearthed sample. The density of soil can be obtained by measuring its volume via the gravity bottle method and making calculations using the relevant formula. A week before the wheat is sown, samples are collected from the ground. The extracted soil sample was cut into appropriate sizes and weighed to obtain the quality of the unearthed sample. Then, the soil samples were placed into a DGF30/7-IA type electric blast (Nanjing experimental instrument factory, Nanjing, China) drying oven at a temperature of 105 °C, baked for 24 h, and weighed.

2.3. Physical Property Parameters of Soil

The pressure subsidence of soil is an important basis for the elastic–plastic deformation and failure of soil [22], which can reflect the ability of soil to withstand pressure. By using the WDW-10 electronic universal testing machine, the soil depth of the settlement and the maximum crushing force can be measured using the uniaxial compression and unconfined compression tests, respectively, to measure the elastic–plastic condition of soil as the index of simulation calibration. The WDW-10 universal testing machine (Shenzhen Suns Technology Stock Co., LTD., Shenzhen, China) is shown in Figure 2.

Experimental Method

During the uniaxial compression test, an appropriate amount of lubricating oil is applied to the inner wall of the plexigroup cylinder with an inner diameter of 50 mm, an outer diameter of 70 mm, and a height of 150 mm to reduce the influence of friction between the inner wall and the sample. To collect samples for the unconfined compressive strength test, the glass cylinder was divided into two semicircles of the same size along the axis and fixed with a banded penstock. The soil sample was placed into a cylinder for testing, and the universal testing machine applied downward pressure, loading at a constant rate of 8 mm/s [12]. The test was repeated four times, as shown in Figure 3.

The soil sample after the uniaxial compression test was taken out and trimmed. The trimmed soil sample was placed on the universal testing machine and pushed downward, with loading at a constant rate of 0.01 mm/s. The test was repeated four times as shown in Figure 4.

2.4. Parameter Calibration of Discrete Element Soil Model

2.4.1. Simulation Parameter Setting

In this study, the EEPA model, which was used for simulation, requires the input of the following parameters: recovery coefficient, static friction coefficient, rolling friction coefficient, constant pull-out force, surface energy, contact/plasticity ratio, slop exp, tensile exp, and tangential stiff multiplier. In order to reduce the influencing factors of the experiment and reduce the difficulty of parameter calibration, the constant pull-out force was not considered in this experiment. The constant pull-out force was defined as 0 and the slop exp was defined as 1. The basic parameters of soil particles were obtained according to the literature [3,12,23], and the specific value was Poisson’s ratio of 0.3.

A cylinder with a height of 150 mm and an inner diameter of 50 mm was built to simulate a glass cylinder. The simulation parameters were set up as shown in

Table 1. Simulation parameters of single-axis closed compression.

Parameter	Value
Recovery coefficient	0.5
Static friction coefficient	0.7
Rolling friction coefficient	0.5
Poisson’s ratio of contact parts	0.3
Shear modulus of contact component	100
Contact component density	7500

.

During the test, the pressure disc was pushed downward at a constant rate of 8 mm/s, and the return rate was the same when the pressure reached 400 N. The change in displacement was recorded from 200 N to 400 N, denoted as the depth of displacement H. The simulation is shown in Figure 5.

After the uniaxial closed compression simulation test, the plate and cylinder were removed, and the sample was stabilized for 1 s to eliminate the kinetic energy generated when the cylinder was moved. The plate moved downward at a constant rate of 0.01 mm/s until the clod broke. Then, the pressure was recorded at its peak, which is denoted as the maximum crushing force F. The simulation is shown in Figure 6.

2.4.2. Plackett–Burman Test

In this study, the constant pulling force was defined as 0, the slop exp as 1, the recovery coefficient as 0.3–0.7, the static friction coefficient as 0.2–1, the rolling friction coefficient as 0.1–0.6, the surface energy as 4–20, the contact/plasticity ratio as 0.3–0.7, the tensile exp as 0.5–3, and the tangential stiff multiplier as 0.3–1.

The parameters that have a significant influence on soil subsidence and maximum pressure values among the parameters to be calibrated can be obtained depending on the range of parameters. Seven parameters needed to be calibrated as influencing factors, and the depression depth and maximum crushing force needed to be calibrated as indexes. The Plackett–Burman test design was carried out using Design-Expert. The test scheme and results are shown in

Table 2. Plackett–Burman test plan and results.

No.	Recovery Coefficient	Static Friction Coefficient	Rolling Friction Coefficient	Surface Energy	Contact/Plasticity Ratio	Tensile Exp	Tangential Stiff Multiplier	Depression Depth	Maximum Crushing Force
1	0.3	0.2	0.1	20	0.2	3	1	2.88	255.5
2	0.3	0.2	0.1	4	0.2	0.5	0.3	3.11	217.0
3	0.3	0.2	0.6	4	0.7	3	0.3	3.12	257.3
4	0.7	1	0.6	4	0.2	0.5	1	3.87	229.5
5	0.3	1	0.1	20	0.7	0.5	1	2.22	266.2
6	0.3	1	0.6	4	0.7	3	1	3.61	252.5
7	0.7	0.2	0.1	4	0.7	0.5	1	3.70	254.5
8	0.3	1	0.6	20	0.2	0.5	0.3	3.81	251.5
9	0.7	0.2	0.6	20	0.7	0.5	0.3	3.21	271.5
10	0.7	1	0.1	20	0.7	3	0.3	3.67	245.5
11	0.7	1	0.1	4	0.2	3	0.3	3.66	223.1
12	0.7	0.2	0.6	20	0.2	3	1	3.70	258.5

.

2.4.3. Central Composite Test

To calibrate the depression depth corresponding to the actual value, Plackett–Burman test factors and a central composite test were used. The test scheme and results are shown in

Table 3. Central composite test scheme and results.

No.	A: Recovery Coefficient	B: Static Friction Coefficient	Depression Depth (mm)
1	0.3	1	3.52
2	0.782843	0.6	3.74
3	0.5	0.6	3.36
4	0.217157	0.6	3.06
5	0.5	0.6	3.57
6	0.5	0.6	3.62
7	0.5	0.6	3.76
8	0.7	1	3.91
9	0.7	0.2	2.67
10	0.5	0.0343146	2.74
11	0.5	1.16569	3.85
12	0.5	0.6	3.72
13	0.3	0.2	3.64

. In the simulation experiment, the middle value calculated using the Plackett–Burman test was used for other non-significant parameters. Specifically, the rolling friction coefficient was 0.5, the surface energy was 10, the plastic deformation ratio was 0.4, the tensile exp was 2, and the tangential stiff multiplier was 0.7.

2.4.4. Box–Behnken Design

To calibrate the maximum destructive power corresponding to the actual value based on Plackett–Burman test factors and the quadratic orthogonal rotation combination test, the Box–Behnken test was carried out. The test scheme and results are shown in

Table 4. Box–Behnken test scheme and results.

No.	C: Surface Energy	D: Contact/Plasticity Ratio	E: Tangential Stiff Multiplier	Maximum Crushing Force (N)
1	12	0.3	0.3	127.61
2	4	0.3	0.65	143.96
3	12	0.5	0.65	174.20
4	4	0.5	1	270.00
5	20	0.3	0.65	256.12
6	12	0.5	0.65	180.40
7	4	0.5	0.3	158.80
8	12	0.3	1	244.81
9	12	0.5	0.65	286.42
10	20	0.7	0.65	416.39
11	12	0.5	0.65	304.38
12	4	0.7	0.65	331.61
13	20	0.5	0.3	267.23
14	20	0.5	1	293.60
15	12	0.7	1	421.60
16	12	0.7	0.3	354.20
17	12	0.5	0.65	294.41

. In the simulation experiment, the recovery coefficient and static friction coefficient are the optimal parameters in the quadratic orthogonal rotation combination test. For other non-significant parameters, the middle value of the Plackett–Burman test was used. Specifically, the rolling friction coefficient was 0.5, the recovery coefficient was 0.322, the static friction coefficient was 0.676, and the tensile exp was 2.

2.4.5. Tangential Stiff Multiplier Optimization Test

The tangential stiffness factor was considered to have a slightly significant contribution to the two indexes of the collapsible depth and maximum crushing force. This was performed to ensure the accuracy of the final data. The single-factor experiment was carried out with the tangential stiff factor as the influence factor, while the depression depth and maximum crushing force were used as indexes. The test scheme and results are shown in

Table 5. Single-factor test scheme and results.

No.	Tangential Stiff Multiplier	Depression Depth (mm)	Maximum Crushing Force (N)
1	0.3	3.28	237.75
2	0.4	3.44	257.25
3	0.5	3.53	236.25
4	0.6	3.41	166.50
5	0.7	3.45	190.50
6	0.8	3.30	201.00
7	0.9	3.54	155.25
8	1	3.61	138.90

.

3. Results

3.1. Analysis of the Intrinsic Parameters of Soil

The density and water content experiments were repeated five times. The results are shown in

Table 6. Results of density and moisture content.

Argument	No.	Result	Average
Density (kg/m3)	1	1629.4	1829.6
	2	1959.2
	3	2152
	4	1370.3
	5	2037.2
Moisture Content (%)	1	22.835	22.84
	2	23.198
	3	22.952
	4	22.262
	5	22.953

. In terms of average values, the density of soil was 1829.6 kg/m³, and the moisture content of soil was 22.84%.

Due to the difference between moisture volatilization and sunshine time caused by soil depth, individual results in the table are either too large or too small. When sowing rotary tillage, the rotary tillage depth can reach 10 cm; therefore, the density and water content are averaged here.

3.2. Analysis of the Physical Property Parameters of Soil

The displacement change curve of uniaxial compression at 400 N is shown in Figure 7, which shows that the growth rate of the compressive force gradually increases. When the pressure is greater than 100 N, the rising rate is relatively stable, and the curve tends to become straight. Therefore, the displacement change, which corresponds to the pressure change from 200 N to 400 N, is denoted as the depression depth H. The final average value is 3.36 mm.

The soil breakage curve is shown in Figure 8. As can be seen from the figure, the curve first gradually increased, rapidly increased to the highest point, slowly decreased, and, finally, the trend gradually leveled off. The highest point of the curve is the maximum pressure when the soil is pressed to complete destruction, which is recorded as the maximum crushing force. The mean value of the maximum crushing force is 206 N.

3.3. Analysis of Parameter Calibration of Discrete Element Soil Model

3.3.1. Plackett–Burman Test and Analysis of Variance

Design-Expert was used to analyze the contribution of compression depth and unconfined pressure obtained using the Plackett–Burman test. The results of the contribution are shown in

Table 7. Contribution degree of depression depth.

Term	Stdized Effect	Sum of Squares	Contribution (%)
Recovery coefficient	1.13173	3.82572	26.9329
Static friction coefficient	1.07517	3.45291	24.3083
Rolling friction coefficient	0.296261	0.262168	1.84565
Surface energy	−0.0597948	0.0106796	0.075184
Contact/plasticity ratio	−0.817339	1.99541	14.0476
Tensile exp	−0.512261	0.783813	5.51801
Tangential stiff multiplier	0.872572	2.27421	16.0104

and

Table 8. Contribution degree of maximum crushing force.

Term	Stdized Effect	Sum of Squares	Contribution (%)
Recovery coefficient	−18.8348	1140.75	0.989614
Static friction coefficient	10.0884	5.80993	4.85739
Rolling friction coefficient	−9.50146	310.083	7.85588
Surface energy	−78.1652	18,018.7	34.7125
Contact/plasticity ratio	185.99	103045	29.7154
Tensile exp	36.3348	4107	0.0601533
Tangential stiff multiplier	64.5015	12,840.1	7.25619

. As can be seen from Table 7, the contributions to the depth of the compressive subsidence (in descending order) are as follows: recovery coefficient, static friction coefficient, tangential stiff multiplier, contact/plasticity ratio, tensile exp, rolling friction coefficient, and surface energy. As can be seen from Table 8, the contributions to the maximum breaking force (in descending order) are as follows: surface energy, contact/plasticity ratio, rolling friction coefficient, tangential stiff multiplier, static friction coefficient, recovery coefficient, and tensile exp. Factors with a contribution degree greater than 15% are considered to be factors with a large contribution degree. It was observed that the recovery coefficient and static friction coefficient significantly contribute to the depth of the compressive subsidence. The surface energy, contact/plasticity ratio, and tangential stiff multiplier have greater contributions to the maximum destructive force.

3.3.2. Central Composite Test and Analysis

The confidence of 0.05 was compared with the p-values of each factor to evaluate the significance of each parameter. As can be seen from

Table 9. Analysis of variance of settlement depth.

Source	Sum of Squares	df	Mean of Squares	F-Value	p-Value
Model	16.81	7	2.4	11.8	0.0076
Recovery coefficient	2.12	1	2.12	10.4	0.0233
Static friction coefficient	5.53	1	5.53	27.15	0.0034
Residual	0.8720	4	0.218
Cor total	17.83	12

, both the recovery coefficient and the static friction coefficient are significant regarding the depth of settlement. With a settlement depth of 3.36 mm as the target, the optimum parameters of the recovery coefficient and static friction coefficient can be obtained. The specific values are as follows: recovery coefficient 0.322, and static friction coefficient 0.676.

Figure 9 shows the influence of the interaction between the recovery coefficient and static friction coefficient on the depth of settlement. Compared with the direction of the static friction coefficient, the surface curve of the recovery coefficient is steeper. The recovery coefficient contour density is slightly higher than the density moving along the static friction coefficient. The effect of the recovery coefficient on the settlement is more significant than that of the static friction coefficient.

According to Table 3, the following regression equation can be obtained by using Design-Expert 12 as follows:

Depth of settlement = 21.587 − 37.179 × A − 39.749 × B + 128.190 × A × B + 19.993 × A² + 25.887 × B² − 71.561 × A² × B − 78.735 × A × B² + 42.199 × A² × B².

A—recovery coefficient.

B—static friction coefficient.

3.3.3. Box–Behnken Design and Regression Analysis

The confidence of 0.05 was compared with the p-values of each factor to evaluate the significance of each parameter.

Table 10. Analysis of variance of maximum crushing force.

Source	Sum of Squares	df	Mean of Squares	F-Value	p-Value
Model	1.062 × 105	9	1.180 × 104	4.58	0.0287
Surface energy	1.353 × 104	1	1.353 × 104	5.25	0.0557
Contact/plasticity ratio	7.522 × 104	1	7.052 × 104	27.36	0.0012
Tangential stiff multiplier	1.299 × 104	1	1.299 × 104	5.04	0.0596
Residual	1.685 × 104	4	4213.55
Cor total	1.242 × 105	16

shows that surface energy, plastic deformation ratio, and tangential stiffness factors have significant effects on the maximum crushing force. The optimal parameters of the surface energy, plastic deformation ratio, and tangential stiffness factor can be obtained using the maximum crushing force of 206 N as the target. The specific values are as follows: surface energy, 17.158; plastic deformation ratio, 0.358; and tangential stiffness factor, 0.447.

Figure 10 shows the influence of the interaction between surface energy, plastic deformation ratio, and tangential stiffness factors on the depth of subsidence. Compared with the direction of the surface energy and tangential stiffness factor, the response surface curve of the plastic deformation ratio is steeper, and the density of the plastic deformation ratio contour line is slightly higher than that moving along the surface energy and tangential stiffness factor. The effect of the plastic deformation ratio on settlement is more significant than that of the static friction coefficient.

According to Table 4, the following regression equation can be obtained by using Design-Expert 12 as follows:

Maximum crushing force = −1.763 + 0.024 × C + 7.037 × D + 5.257 × E − 0.783 × C × D + 0.138 × C × F − 4.429 × D × E + 0.026 × C² + 3.494 × D² + 0.424 × E² − 0.037 × C² × D − 0.014 × C² × E + 1.636 × C × D².

C—surface energy.

D—contact/plasticity ratio.

E—tangential stiff multiplier.

3.3.4. Analysis of Tangential Stiffness Factor Optimization Test Results

Taking the actual values of the compression depth (3.36 mm) and maximum crushing force (206 N) as target lines, a comparison of single-factor experimental results is shown in Figure 11. The errors between the depth of compression and maximum crushing force and the actual values are calculated as shown in Figure 12. It can be intuitively seen from the figure that the error is the smallest at 0.8, and thus the tangential stiffness factor is 0.8.

According to the above test and result analysis, the optimal combination of parameters to be calibrated can be obtained with the following parameters: recovery coefficient, 0.322; static friction coefficient, 0.676; rolling friction coefficient, 0.5; surface energy, 17.158; contact/plasticity ratio, 0.358; tensile exp, 2; tangential stiff multiplier, 0.8.

4. Verification of Calibration Results

According to the calibration results of the depth of settlement and maximum crushing force parameters, a discrete element simulation model was established using EDEM(2020, DEM-Solutions, Edinburgh, UK) software, and uniaxial compression and unconfined compression tests were carried out to verify the reliability of the calibration results. The actual and simulated values of the uniaxial compression test curve are shown in Figure 9. The actual and simulated values of the unconfined compressive test are shown in Figure 10.

OriginPro( 2021, OriginLab, Northampton, MA, USA) software was used to analyze the two curves obtained from the uniaxial compression experiment as well as the unconfined compression test simulation and actual test, respectively. Figure 13 shows the curve comparison between the simulation and actual test results of the uniaxial compression test. The analysis results show that the error between the actual value, where the simulation value is 2.782%, and the misalignment coefficient R² is 93.509%. Figure 14 shows the curve comparison between the simulation results of the unconfined compressive test and the actual test results. The analysis results show that the error between the actual value and the simulated value is 3.674%, and the misfit coefficient R² is 94.2%. The analysis results of these two graphs show that the two curves of the actual value and the simulation value are highly similar. Therefore, the discrete element soil model after parameter calibration can be used for the simulation optimization of the subsequent mechanism.

5. Discussion

In recent years, it has become increasingly common to use simulation software to simulate field seeding and growth in agriculture and related fields. Simulations are not subjected to the restrictions of weather, temperature, season, and other changes that affect field testing. The most difficult challenge is selecting the correct contact model of soil particles and the appropriate calibration method for the model’s parameters. Wang Xianliang selected the EEPA model and the plate subsidence as the index, screened four factors (recovery coefficient, static friction coefficient, surface energy, contact/plasticity ratio) using the Plackett–Burman test, and then obtained the optimal solution using bbd. Janda [19] presents a numerical study on the DEM modeling of cohesive solids using a visco-elastoplastic frictional adhesive contact model. The capabilities of the contact model to capture the mechanical macroscopic behavior of cohesive materials were investigated by means of cone penetration and unconfined compression simulations. Qishuo Ding [24] selected the Hertz–Mindlin model with a bonding soil particle contact model, considered field-measured tillage resistance as the standard value for verifying critical stress, and established a subsoiling shovel–soil discrete element tillage model suitable for the test process. Based on the discrete element method, Junwei Li [17] used the Hertz–Mindlin model alongside the JKR cohesion model to conduct a discrete element simulation of two kinds of viscous black soils with different moisture contents. Through the simulation test of the accumulation angle, the simulation parameters of the contact model between soil particles under the two values of moisture content were calibrated. The parameters of the two soil contact models with water content were obtained using a response surface optimization method. Xie Fangping measured the compressive strength of soil under pressure using an unconfined compressive strength test and obtained the results that best fit the target value using pb and bbd tests. Wanli Pengcheng designed a stacking angle test and a compression test of the single axle seal to obtain the most suitable results.

When comparing the experimental results of different scholars, it is found that the depression depth curve obtained via the uniaxial compression test has the same trend as the pressure curve obtained by Xianliang Wang. The pressure slowly increases at first before rapidly increasing and finally showing a linear sloping trend. The trend of the maximum crushing force curve is similar to that determined by Fangping Xie and Janda. The curve first rapidly rises to its highest point, rapidly declines for a while, and finally slowly declines to a straight horizontal line. However, the maximum crushing force is slightly larger than the force determined by Fangping Xie and Janda since the soil sampled in this experiment has a larger moisture content and is slightly more viscous than the experimental soil found by other experts and scholars.

In this study, the parameters of rice–stubble wheat soil in the Yangtze River Basin of China were calibrated with the index of the collapsing depth and maximum crushing force. This study promotes the development of a soil model for wheat grown after rice in Jiangsu, as well as promoting the optimization of sowing machine components for growing wheat after rice. This paper provides the methods and data necessary for conducting simulations on soil used for wheat grown after rice in the Yangtze River Basin of China. However, there are still some major shortcomings, such as the following:

(1) The static friction and collision recovery coefficients of parameters to be simulated can be measured using a test. However, the two parameters used in this study refer to previously investigated research. For the accuracy of the calibration value, accurate values should be measured first to reduce the number of calibration parameters.

(2) We studied the soil for wheat grown after rice in the Jiangsu Province. Different soils in different areas of the Yangtze River Basin have different moisture contents and viscosity levels, and models should be constructed according to these specific conditions.

6. Conclusions

In this study, EDEM (2020, DEM-Solutions, Edinburgh, UK) software was used to design and analyze the EEPA model via the Plackett–Burman test, quadratic orthogonal rotation combination test, Box–Behnken test, and single-factor test schemes using Design-Expert. The soil simulation parameters for wheat grown after rice in Jiangsu were calibrated. Comparing the results of both the parameter calibration test and analysis, the main conclusions are as follows:

(1) The Plackett–Burman test results show that the recovery coefficient and static friction coefficient of the model’s parameters have a large influence on the depth of compression. The surface energy, contact/plasticity ratio, and tangential stiff multiplier have a large influence on the maximum destructive force.

(2) The results of the quadratic orthogonal rotation combination test, Box–Behnken test, and single-factor test are as follows. The recovery coefficient corresponding to the actual value is 0.322, the static friction coefficient is 0.676, the rolling friction coefficient is 0.5, the surface energy is 17.158, the contact/plasticity ratio is 0.358, the tensile exp is 2, and the tangential stiff multiplier is 0.8.

(3) The calibration results show that the error between the actual value and the simulation value of the uniaxial compression and unconfined compression tests is less than 5%, and the misfit coefficient R² is greater than 90%. It can be seen that the model and the results obtained can be applied to the optimization simulation of the subsequent mechanism with a high accuracy.