Calibration Parameter of Soil Discrete Element Based on Area Difference Method

Abstract

Existing calibration methods for contact parameters of soil discrete elements are primarily based on a single physical property of soil (e.g., natural accumulation angle and shear characteristics). To solve this limitation, a calibration method based on the area difference method was proposed in accordance with the actual operation effect of ridge and furrow opener. The typical clay loam in the middle and lower reaches of the Yangtze River (moisture content of 20.66% ± 1%) was used as the research object. The Hertz–Mindlin with JKR Cohesion contact model was selected, and the static friction factor and rolling friction factor between soil and 65Mn steel were determined by inclined plane mechanics method. Soil tank ridge and furrow opening experiments and quadratic orthogonal rotation combined tests were carried out. MATLAB software was used to carry out polynomial fitting and the constant integral solution to the ridge and the cross-sectional contour data of the furrow obtained from the measured and simulated ridge and furrow opening experiments to establish the regression model of area difference. The optimization function in Design-Expert software was adopted. The optimal parameter combination was obtained by optimizing the regression model with the minimum area difference as the target. Based on the above parameter combination, the relative error of cross-sectional area of the furrow that was obtained by simulation test and soil tank test was 9.15%. Finally, the comparison of the traction resistance of the simulation and the soil tank test at the stable working stage indicated that the relative error between the simulated value and the measured value was 10.4%, which is less than the acceptable error of 15% in general statistics and within a reasonable range. The results showed that the particle movement and mechanical properties of the soil simulation model were basically consistent with the actual soil after parameter optimization and calibration. In this study, a precise discrete element simulation model of the typical clay loam in the middle and lower reaches of the Yangtze River based on the area difference method was constructed, which lays a theoretical basis for later exploration of the interaction mechanism between the typical clay loam and the soil contact components in the middle and lower reaches of the Yangtze River.

1. Introduction

Research on the interaction law between soil and soil contact components is the basis of design [1]. The discrete element method is based on the dynamic relaxation method, Newton’s second law and time-step iteration to solve the motion velocity and displacement of each particle; it is capable of analyzing the force and motion law of soil particles under the action of soil contact components from a microscopic perspective [2]. The discrete element method has served as a general method to study the dynamics of loose particles and has been extensively used to analyze the interaction law between soil and soil contact components [3,4,5,6].

The calibration of the simulation parameters is the key to constructing a reliable discrete element simulation model, which has been explored by experts and scholars worldwide based on the combination of actual measurements and simulation tests to better calibrate discrete element simulation parameters. Aikins et al. [7] used a combination of measured tests and discrete element simulation to calibrate contact parameters for high moisture content clay loam through stacking experiments. They verified the results by comparing the effect of narrow point opener geometry on soil disturbance and cutting forces in simulated and measured tests. Wang et al. [8] substituted the Linear Cohesion contact model into the Hysteretic Spring contact model to simulate soil cohesion and calibrated the discrete element simulation parameters for loose clay soils. Given the plastic deformation between soil particles, Ucgul et al. [9] used the Hysteretic Spring contact model to calibrate the friction factor between soil particles by stacking angle and the direct shear test. Coetzee. [10] determined the particle stiffness using a lateral limit compression test, and the interparticle friction factor was optimally calibrated through a direct shear test and a stacking angle test. Xiang et al. [11] optimally calibrated discrete element simulation parameters for clay loam in southern China through stacking experiments, and the accuracy of the simulation parameters was verified using cavity formation tests. Wang et al. [12] used the Edinburgh Elasto-Plastic Cohesion Model based on the agent model approach, and the particle radius, static and rolling friction factors of the soil were optimally calibrated through stacking experiments and direct shear test; the effectiveness of the method was verified using tire–soil compaction tests. Li et al. [13] optimally calibrated the contact parameters between the clayey black soil particles and the contact parameters of the clayey black soil with 65Mn steel and PTFE plates, respectively. Sun et al. [14] calibrated the simulated contact parameters for clay loam on typical slopes of the Chinese Loess Plateau, and the reliability of the simulated parameters was verified by slope rotational plowing tests. Song et al. [15] calibrated the simulation parameters between soil particles and between soil and soil contact components in a post-tillage cotton field using stacking experiments and slip tests, and the optimally calibrated contact parameters were verified by comparing the working resistance between field and simulation tests during furrow opening and soil covering processes.

The above methods of calibration between soil particles’ contact parameters are primarily based on some single physical characteristics of soil (e.g., natural accumulation angle and shear characteristics). In this study, a method of calibration for soil discrete element simulation parameters based on actual operational effects is proposed. For the typical clay loam in the middle and lower reaches of the Yangtze River, the Hertz–Mindlin with JKR Cohesion contact model and Design-Expert software were used in this study. According to the MATLAB software (Version No. 2016a, MathWorks Inc., Neddick, MA, USA), polynomial fitting and definite integral solving were carried out on cross-sectional contour data of the furrow that were obtained by the moldboard-style ridge and furrow opener. The resulting area difference was the response value. A Box–Behnken test was designed to optimize and calibrate the coefficient of recovery, static friction, rolling friction, and soil JKR surface energy. According to the calibrated parameters, the discrete element simulation model of the typical clay loam in the middle and lower reaches of the Yangtze River was constructed. Finally, the reliability of the calibration parameters was verified by simulation test and traction resistance verification test, which lays a theoretical basis for the subsequent simulation and optimization design of typical clay loam and earth-touching components in the middle and lower reaches of the Yangtze River.

2. Materials and Methods

2.1. Determination of Cross-Sectional Contour Data of the Furrow

2.1.1. Test Conditions

The experiment was conducted in a digital soil tank at the Agricultural Mechanization Engineering Training Center of Hunan Agricultural University, and the soil in the tank originated from the experimental field of the Cultivation Garden of Hunan Agricultural University. First, the particle size distribution and texture composition of the soil tank were analyzed using the screening method. The screen was used for screening, and then the electronic balance with the precision of 0.01 g and the range of 0~300 g was used to weigh the soil. The test was repeated three times, and the particle size distribution and percentage content of the soil samples were obtained: d > 0.9 mm (14.6%), d = 0.9~0.77 mm (3.2%), d = 0.77~0.5 mm (16.7%), d = 0.5~0.4 mm (12.7%), d < 0.4 mm (52.8%). According to the Chinese soil geological classification standard [16], the soil type in the soil tank is typical clay loam in the middle and lower reaches of the Yangtze River. Considering that the soil tank is located indoors, soil samples with soil depths of 0~200 mm were obtained in the experimental field of the Cultivation Garden of Hunan Agricultural University by the five-point sampling method, and the soil moisture content of 20.66% was examined using a drying experiment, and the respective sampling point was repeatedly examined three times to obtain the average value, consistent with the actual moisture content of South China production. To ensure that the soil moisture content of each part of the soil tank was consistent with the soil moisture content measured in the test field of Hunan Agricultural University’s Cultivation Garden, the soil moisture content within the test range of the soil tank was determined by the five-point sampling method based on a hand-held soil moisture content rapid tester. If the soil moisture content of the soil tank was examined to be lower than 20.66% ± 2%, the soil was poured with an appropriate amount of water using a sprinkler and stood for 3~5 h. Subsequently, the soil was rototilled at a depth of 0~25 cm using a rototiller, such that the upper and lower layers of soil were mixed evenly, and the soil water content was examined. If the measured soil moisture content of the soil tank was higher than 20.66% ± 2%, the rototiller was adopted to rototill the soil at 1 h intervals and measure the soil moisture content, and the above steps were repeated. Moreover, the soil moisture content of the test tank was adjusted to stabilize at 20.66% ± 2%. The soil compactness of 0~476.5 KPa at depth 0~200 mm was examined using a TYD-2 digital display soil hardness meter (Ai&Ce Inc., Zhejiang, China), and the soil density of 2.57 g/cm³ was examined through a specific gravity bottle experiment, thus laying the basis for the subsequent discrete element soil simulation model establishment.

2.1.2. Ridge and Furrow Opening Test

The TCC electric four-wheel drive soil tank test-bed is presented in Figure 1a, with the moldboard-style ridge and furrow opener (Figure 1b) mounted on the rear beam of the sextant force measurement frame. The real-time data collector of the sextant force measurement frame was connected to the computer console of the soil tank test-bed through a data cable, and the sextant force measurement frame was hooked up to the soil tank test-bed through a translatable and liftable three-point suspension device.

The control platform of the soil tank test-bed was opened, and the trench was dug out under the ridge and furrow opener, whose length and width were slightly greater than the opener and whose depth was consistent with the theoretical furrow-type contour (Figure 2). The intersection line of the upper and lower wings of the opener was adjusted to be flush with the ground through the control platform. The soil tank test-bed measurement and control software were opened to store the traction resistance of the opener, the running speed of the soil tank test-bed was set to 0.72 m/s (the speed was obtained from the speed of the Dongfang Hong-MG600 tractor at full load, i.e., low-speed first gear 2.6 km/h), and the soil tank test-bed was run to the end of 30 m in front and stopped. To ensure the accuracy of the test data, the acceleration and deceleration sections of the soil tank test-bed were reduced by 7 m, the cross-sectional contour parameters of the middle 16 m long furrow were examined using a furrow contour measuring device (Figure 3) and the data were measured at 1 m intervals for a total of 16 groups.

2.1.3. Data Processing and Analysis

The data of 16 sets of cross-sectional contour data of the furrows measured in the experiment were averaged, and the two-dimensional data coordinate points obtained from the test data were geometric translational transformations to generate a point line diagram of the cross-sectional contour parameters of the furrows (Figure 4), and the soil simulation parameters using the area difference method were calibrated with the results obtained from this test as the target value. The data of traction resistance in the forward direction of the ridge and furrow opener measured in this test were prepared for the subsequent traction resistance verification test, which will not be analyzed in depth here.

2.2. Discrete Element Simulation Modeling

2.2.1. Contact Model Selection

A typical clay loam in the middle and lower reaches of the Yangtze River was used in this study. The soil in this area has high moisture content and strong cohesion and elastoplasticity between soil particles. The Hertz–Mindlin (no slip) contact model is accurate and efficient in the calculation of forces, whereas it does not consider the effect of interparticle bonding forces on the particle motion law [17]. The Hertz–Mindlin with Bonding contact model is suitable for bonding particles, whereas the particles are only adhered together by limited size bonding bonds. When the bond breaks under external forces, the separated particles will exist independently of each other [18]. Neither of the above two contact models can accurately simulate the dynamic motion law of viscous soil particles under the action of soil contact components. The Hertz–Mindlin with JKR Cohesion contact model is designed to simulate cohesive wet particles and can fully consider the effect of the cohesive force between soil particles on the particle motion law. Thus, this contact model was selected to perform simulation to ensure the accuracy of the calculation results.

2.2.2. Determination of Soil–Steel Friction Factor

The material of the ridge and furrow opener used in this study was 65Mn steel. The friction factor between the soil and the soil contact components affects the working effect and resistance characteristics of the contact parts. To ensure the accuracy and reliability of the calibration parameters, the static friction factor and rolling friction factor between the soil and 65Mn steel were determined using the inclined plane mechanics method [19,20,21].

The 65Mn steel plate was pretreated before the test, the rust and residual debris on the surface of the plate was polished until smooth and flat, so as not to affect the test results. An inclined plane test bench was built (Figure 5a). The center position of the thin rope and SANS-CMT6104 universal mechanical testing machine (SANS Inc., Shenzhen, Guangdong, China), through the mechanics of the probe connected to the two ends of the thin rope symmetrically tied any narrow side to the 300 mm × 500 mm 65Mn steel plate. VICTOR digital display inclinometer was placed along the length of the steel plate, the vertical height of the universal mechanical testing machine mechanics probe was adjusted until the digital display inclinometer showed 0 degrees to ensure the level of the steel plate. Soils with 20.66% (±1%) moisture content were shaped.

(1)

Determination of Soil–Steel Static Friction Factor.

The soil with 20.66% (±1%) moisture content was shaped into a cylindrical soil block; the principle diagram of the inclined plane mechanics method is presented in Figure 5b, and the force analysis of the soil block was conducted; air resistance, airflow changes and other factors were ignored.

$$\{\begin{array}{l}{f}_1=m_{1}g\cos \alpha \\ N_1=m_{1}g\sin \alpha \\ f_1=\mu f{N}_1\\ \mu =\tan \alpha \end{array}$$

(1)

where f₁ is the frictional force on the soil block, N; N₁ is the support force on the soil block, N; α is the angle between the 65Mn steel plate and horizontal surface, (°); μ is the static friction factor of the soil block and 65Mn steel plate.

The soil block was placed in the center of the steel plate width direction, and the mechanical probe of the mechanical testing machine was raised at the speed of 10 mm/min. Once the clod began to slide on it, the inclination of the incline was stopped from increasing. The VICTOR digital display inclinometer was used to record the dip angle of steel plate at this time, and three repeated tests were conducted to obtain the mean value. The test results are shown in

Table 1. Soil–steel static friction factor measurement results.

Serial Number	65Mn Steel Plate Inclination Angle/(°)	Static Friction Factor	Mean Value
1	34.2°	0.68	0.71
2	35.3°	0.71
3	36.4°	0.74

.

(2)

Determination of Soil–Steel Rolling Friction Factor.

First, the soil with a moisture content of 20.66% (±1%) was made into a soil ball with a radius of 10 mm; the soil–steel rolling friction factor was determined the same use of inclined plane mechanics method, and the force analysis was conducted under the neglect of air resistance, airflow changes and other factors (Figure 6).

$$\{\begin{array}{l}{f}_2=m_{2}g\cos \beta \\ N_2=m_{2}g\sin \beta \\ f_2=\mu f{N}_2\\ k=\tan \beta \end{array}$$

(2)

where f₂ is the frictional force on the soil ball, N; N₂ is the support force on the soil ball, N; β is the angle between 65Mn steel plate and horizontal surface, (°); k is the rolling friction factor of the soil ball and 65Mn steel plate.

Moreover, the soil ball was placed in the center of the width of the steel plate, and the mechanics probe of the mechanics tester was lifted at a speed of 10 mm/min. When the soil ball was about to roll on the steel plate, the mechanics probe was stopped from being lifted. The digital display inclinometer was adopted to record the inclination angle of the inclined plane at this time. The test was repeated three times to obtain the average value, and the test results are listed in

Table 2. Soil–steel rolling friction factor measurement results.

Serial Number	65Mn Steel Plate Inclination Angle/(°)	Rolling Friction Factor	Mean Value
1	5.6°	0.098	0.097
2	5.2°	0.091
3	5.8°	0.102

.

2.2.3. Calibration Parameters Selection

In the discrete element simulation modeling process, some of the soil intrinsic parameters can be examined directly through physical experiments (e.g., soil moisture content and soil particle density), whereas there are also some soil contact parameters that are more difficult to measure accurately in practice (e.g., soil–soil static friction factor, soil–soil rolling friction factor, soil collision recovery coefficient, and soil JKR surface energy). Thus, the above contact parameters should be calibrated through simulation optimization tests. Before the simulation optimization tests, the value range of the respective calibration parameter should be determined. Over the past few years, the application of discrete element simulation in the field of agricultural engineering has become increasingly advanced. Notably, the research and application of the mechanism of interaction between soil contact components and soil has achieved more results. In this study, the value ranges of each calibration parameter were determined with full reference to the existing research results [11,12,13,22,23,24,25,26]: recovery coefficient: 0.3~0.6, static friction factor: 0.2~0.6, rolling friction factor: 0.05~0.3 and JKR surface energy: 5 J/m²~9 J/m².

2.2.4. Establishment of the Simulation Test Model

In actual soils, soil particles are generally classified in accordance with size as stone, sand–gravel, powder and clay particles. The range of particle sizes is respectively: 2~10 mm, 0.05~1 mm, 0.005~0.1 mm and 0.0001~0.01 mm [27], with varying particle shapes. It is unlikely to model soil particles with a 1:1 ratio based on their actual triaxial size and shape. Existing research has suggested that in the discrete element simulation of soil contact components and soil interaction process, the smaller the radius of the soil particle model, the more accurate the simulation calculation results will be, whereas the simulation time will also increase. When the soil particle radius is 10 mm, the calculation results can still accurately indicate the actual operation of soil contact components during the force and operation effect [28]. In this study, to ensure the accuracy and efficiency of the simulation test simultaneously, a spherical soil particle model with a radius of 8 mm was built, while a soil particle model with 0.9~1.1 times the basic particle radius was randomly generated by setting in the particle plant. The measurement of soil moisture content and density parameters has been completed in the previous research, and the soil Poisson’s ratio ν = 0.38 and shear modulus G = 1 MP were selected based on the physical properties of typical clay loam in the middle and lower reaches of the Yangtze River and a review of the relevant literature [29]. The ridge and furrow opener was modeled in three dimensions using SolidWorks software (Version No. 2016, Dassault Systemes Inc., Concord, Massachusetts, USA) at a ratio of 1:1, saved in IGS format, and substituted into EDEM software (Version No. 2018, DEMSL Inc., Edinburgh, Scotland, UK). Hertz–Mindlin (no slip) was selected as the contact model between the ridge and furrow opener and the soil. The intrinsic physical parameters of 65Mn steel were reviewed in the literature [30]: Poisson’s ratio 0.30, density 7861 kg/m³, shear modulus 7.9 × 10¹⁰ Pa; soil–65Mn steel recovery coefficient 0.3 [11].

A simulated soil tank model of 3000 mm × 1200 mm × 300 mm (L × W × H) was built. Moreover, the advancing speed of the ridge and furrow opener was set to 0.72 m/s, and the depth of soil entry was set to 200 mm. The Rayleigh time step was obtained as 20%, the simulation time step was 3.1 + 05 es and the simulation time reached 16 s. The cell grid size was set to 2 times of the minimum particle radius, and the process of simulating the ridge and furrow opening operation is presented in Figure 7.

3. Results and Discussion

3.1. Design of Calibration Experiment

To accurately evaluate the difference between the furrow opening effect of the soil tank ridge and furrow opening experiments and the simulated ridge and furrow opening experiments, based on the cross-sectional area difference of the furrows obtained from the soil tank ridge and furrow opening experiments and the simulated ridge and furrow opening experiments, polynomial fitting was performed on the cross-sectional contour data of the furrow obtained from Section 2.1.3 soil tank ridge and furrow opening experiments using MATLAB2016a [31] software. The polynomial obtained is expressed in Equation 3, and the image of the corresponding furrow-type contour function of F(x) is presented in Figure 8.

$$\begin{array}{c}F\left(x\right)=-5.04\times {10}^{-19}\cdot x^9-1.96\times {10}^{-17}\cdot x^8\\ +5.65\times {10}^{-14}\cdot x^7+2.12\times {10}^{-12}\cdot x^6-2.05\times {10}^{-9}\cdot x^5\\ -9.65\times {10}^{-8}\cdot x^4+2.97\times {10}^{-5}\cdot x^3+4.41\times {10}^{-3}\cdot x^2\\ -0.27x+7.28\end{array}$$

(3)

After the completion of the simulated ridge and furrow opening experiments, an area with a length of 1 m from the central segment of the simulated soil tank was used, and a truncation analysis of the simulated soil tank was conducted along the direction perpendicular to the advancing direction of the ridge and furrow opener (Y-axis direction) using the Clipping module of EDEM software, with a slice thickness of 50 mm and a slice interval of 100 mm, and a total of six slices were obtained. The simulated soil particles were obtained in each slice that can represent the furrow-type contour. As shown in Figure 9a, the two-dimensional scatter plot of the X and Z coordinate data of the particles was generated. Subsequently, the data coordinate points of the particles were transformed by geometric translation to ensure that the graph was essentially axis symmetric at x = 0 to facilitate the subsequent area difference calculation, and the six sets of coordinate data from the above preliminary processing were imported into MATLAB software for polynomial fitting to obtain the polynomial Φ(x) for the cross-sectional furrow-type contour of the furrow. Six polynomials were obtained after the fitting, which will not be listed one by one in this study. Any one of these polynomial fit function images is shown in Figure 9b.

After the function equations of F(x) and Φ(x) were obtained, separately, the area of the enclosed area between F(x) and Φ(x), as shown in Figure 10, was solved by definite integration using MATLAB software, and the area difference of the six slices was calculated separately by obtaining the mean value:

$$\Delta S=\frac{1}{n}{\displaystyle \sum _{i=1}^n\Delta S_i}=\frac{1}{n}{{\displaystyle \sum _{i=1}^n{\displaystyle {\int }_a^b\left[{\Phi }_i(x)-F_i(x)\right]}}}_{i}dx$$

(4)

In other words, the cross-sectional area difference of the furrows was obtained from the soil tank ridge and furrow opening experiments and the simulated ridge and furrow opening experiments (denoted by ΔS).

In this experiment, a four-factor, three-level quadratic orthogonal rotational combination test was conducted with soil-to-soil recovery coefficient, static friction factor, rolling friction factor and soil JKR surface energy as test factors and area difference ΔS as evaluation index, and the test factors and levels are listed in

Table 3. Calibration of test factors and levels.

Experimental Factors	Low Level	Intermediate Level	High Level
Soil–soil recovery coefficient	0.3	0.45	0.6
Soil–soil static friction factor	0.2	0.4	0.6
Soil–soil rolling friction factor	0.05	0.175	0.3
Soil JKR surface energy (J/${\mathrm{m}}^2$)	5	7	9

.

3.2. Simulation Results and Analysis

A three-level quadratic orthogonal rotational combination test was designed in accordance with the Box–Behnken test principle, and the test design scheme and results are listed in

Table 4. Experimental design scheme and results.

Serial Number	Soil–Soil Recovery Coefficient	Soil–Soil Static Friction Factor	Soil–Soil Rolling Friction Factor	Soil JKR Surface Energy/(J/m−2)	Area Difference/(dm2)
1	0.30	0.20	0.17	7.00	2.00
2	0.60	0.20	0.17	7.00	0.56
3	0.30	0.60	0.17	7.00	2.44
4	0.60	0.60	0.17	7.00	1.17
5	0.45	0.40	0.05	5.00	0.89
6	0.45	0.40	0.30	5.00	1.72
7	0.45	0.40	0.05	9.00	0.35
8	0.45	0.40	0.30	9.00	0.87
9	0.30	0.40	0.17	5.00	2.12
10	0.60	0.40	0.17	5.00	1.07
11	0.30	0.40	0.17	9.00	2.13
12	0.60	0.40	0.17	9.00	0.80
13	0.45	0.20	0.05	7.00	0.20
14	0.45	0.60	0.05	7.00	0.70
15	0.45	0.20	0.30	7.00	0.99
16	0.45	0.60	0.30	7.00	1.18
17	0.30	0.40	0.05	7.00	1.78
18	0.60	0.40	0.05	7.00	0.71
19	0.30	0.40	0.30	7.00	2.07
20	0.60	0.40	0.30	7.00	0.64
21	0.45	0.20	0.17	5.00	0.40
22	0.45	0.60	0.17	5.00	1.44
23	0.45	0.20	0.17	9.00	0.58
24	0.45	0.60	0.17	9.00	0.60
25	0.45	0.40	0.17	7.00	0.64
26	0.45	0.40	0.17	7.00	1.03
27	0.45	0.40	0.17	7.00	0.75
28	0.45	0.40	0.17	7.00	0.80
29	0.45	0.40	0.17	7.00	0.87

.

3.2.1. Regression Modeling and Significance Analysis

The quadratic regression analysis of the test results was conducted using Design-Expert 11 software. In accordance with $\Delta {\mathrm{S}}_{\max }/\Delta {\mathrm{S}}_{\min }$ = 11.9 > 10, the square root transformation of S was performed, such that $\Delta {\mathrm{S}}_{\max }/\Delta {\mathrm{S}}_{\min }$ = 3.4. The multiple regression equation of $\sqrt{\Delta \mathrm{S}}$ was established (as expressed in Equation (5)) and tested for significance.

$$\begin{array}{c}{Y}_1=0.90-0.27A+0.13B+0.13C\\ -0.097D+0.048AB-0.072BC\\ -0.14BD+0.28A^2-0.032B^2\end{array}$$

(5)

where A denotes the soil-to-soil recovery coefficient; B is the soil-to-soil static friction factor; C is the soil-to-soil rolling friction factor; D is the soil JKR surface energy; Y₁ represents the $\sqrt{\Delta \mathrm{S}}$.

The results of the analysis of variance of the model are listed in

Table 5. Analysis of variance for $Y_1$ regression model.

Source of Variance	Sum of Squares	Degree of Freedom	Mean Square	F Value	p Value	Significance
Model	2.10	9	0.23	19.97	<0.0001	**
A	0.88	1	0.88	75.60	<0.0001	**
B	0.20	1	0.20	16.92	0.0006	**
C	0.21	1	0.21	17.58	0.0005	**
D	0.11	1	0.11	9.67	0.0058	**
BD	0.079	1	0.079	6.76	0.0176	*
A2	0.56	1	0.56	47.59	<0.0001	**
Residual	0.22	19	0.012
Lack of fit	0.20	15	0.013	2.11	0.2459
Pure error	0.025	4	6.23 × 10-3
Summation	2.32	28
R2 = 0.9044     adj-R2 = 0.8591     C.V. = 10.77%     Adeq Precisior = 16.537

Note: p indicates the probability level criterion to test whether it is significant, p < 0.05 denotes that the significant model term, expressed by “*”, p < 0.01 indicates that the model term is highly significant, “**” represents that p > 0.1 indicates the model term that is not significant.

. The p value of the equation model was <0.0001, which is highly significant (p < 0.01); the coefficient of determination R² = 0.9044 and the corrected coefficient of determination adj-R² = 0.8591, both of which were close to 1; the coefficient of variation C.V. = 10.77% and the precision Adeq Precisior = 16.537; indicating that the model exhibited high reliability and good precision, and the cross-sectional area difference of the furrows obtained from the soil tank ridge and furrow opening experiments and the simulated ridge and furrow opening experiments can be predicted according to this model.

The effects of the one degree term A, B, C and D and the secondary term A² of the model were highly significant and the interaction term BD was significant within the reference values of the above test factors. According to the magnitude of the regression coefficients of the respective factor in the Y₁ regression model, the order of influence of the test factors on the area difference was obtained: A > C > B > D.

3.2.2. The Influence Law of Interaction Factors on Area Difference

The influence law of recovery coefficient A, static friction factor B, rolling friction factor C and JKR surface energy D interaction on the area difference was analyzed in accordance with the results of analysis of variance of the regression model. Moreover, the response surface was plotted using Design-Expert software (Version No. 11, Stat-Ease Inc., Minneapolis, MN, USA), and the response surface of the effect of interaction factors on the area difference ΔS is presented in Figure 11.

As depicted in Figure 11a, when the static friction factor B was certain, the area difference first decreased and then increased with the increase in the recovery coefficient A, and there was an optimal value of the recovery coefficient A in its range of values. When the recovery coefficient A was certain, the area difference decreased with the decrease in the static friction factor B. Moreover, at a certain recovery coefficient A, the static friction factor B was positively correlated with the area difference. As depicted in Figure 11b, simultaneously reducing the static friction factor B and rolling friction factor C was conducive to reducing the area difference, such that the static friction factor B and rolling friction factor C were positively correlated with the area difference. As depicted in Figure 11c, decreasing the static friction factor B and increasing the JKR surface energy D contributed to the reduction of the area difference, such that the static friction factor B was positively correlated with the area difference, and the JKR surface energy D was negatively correlated with the area difference.

3.2.3. Optimal Parameter Groups and Simulation Test Verification

Optimization analysis was performed by the optimization function in Design-Expert software, and the regression model was searched for the minimum value of area difference as the target. The optimal combination of parameters was determined as soil-to-soil recovery coefficient of 0.5, static friction factor of 0.2, rolling friction factor of 0.05 and soil JKR surface energy of 5 J/m².

The optimal combination of parameters was used for simulation test verification, and the ridge and furrow opener applied and other parameters were kept constant. The test was repeated three times, and the area difference was determined as 0.29, 0.29 and 0.32, with a mean value of 0.30 and a standard deviation of 0.02. The relative error of the measured cross-sectional area of the furrow with the previous soil tank test was 9.15%, less than the general statistically acceptable error of 15% and within a reasonable range. The results suggest that the model based on the area difference method can accurately reflect the movement pattern of soil particles.

3.3. Optimal Parameter Groups and Simulation Test Verification

The mechanical properties of the simulated soil model were verified by comparing the traction resistance obtained from the simulated ridge and furrow opening operation with the actual traction resistance obtained from the previous soil tank test under the optimal combination of parameters. A certain degree of compaction was applied to the soil tank model at the early stage of the simulation test to achieve a soil compactness of 0~476.5(±2%) KPa at a depth of 0~200 mm. Moreover, the forward speed of the opener was set to 0.72 m/s with a soil entry depth of 200 mm.

The working process was divided into three stages when the moldboard-style ridge and furrow opener were employed for ridge and furrow opening operation. To be specific, the soil entry stage, the stable working stage and the soil exit stage were set as the first, second and final stages. The traction resistance at the stable working stage was stable at a certain value and fluctuated slightly up and down. The traction resistance curve at the stable working stage was analyzed, as presented in Figure 12. The traction resistance at the stable working stage showed irregular up and down fluctuation with time; the average value of traction resistance at the stable working stage was 1303.07 N and 1167.30 N in the soil tank test and simulation test, respectively, the traction resistance calculated by the simulation soil model was smaller than the traction resistance of actual test. The reason for the above result is that crop stubble and other substances will inevitably exist in the soil during the actual test, thus resulting in an increase in the traction resistance of the opener.

In general, the basic trend of the test curve was consistent with that of the simulation curve, and the relative error of traction resistance between the stable stage of the simulated ridge and furrow opening experiments and the stable phase of the soil tank ridge and furrow opening experiments was 10.4%, less than the general statistically acceptable error of 15%, and the error fell into a reasonable range, thus suggesting that the mechanical properties of the simulated soil model were basically consistent with the actual soil. The accuracy of the discrete element simulation parameters based on the area difference method was verified.

4. Conclusions

In this study, according to the actual furrow opening effect, an optimal calibration method for soil discrete element simulation parameters based on the area difference method was proposed, and the accuracy of the area difference method was verified by the precision of the discrete element simulation model of a typical clay loam soil in the middle and lower reaches of the Yangtze River constructed by the Hertz–Mindlin with JKR Cohesion contact model. Firstly, the soil–steel contact parameters were determined using the inclined plane mechanics method; the area difference regression model was built by the soil tank ridge and furrow opening test and the quadratic orthogonal rotation combination test, and the significance analysis showed that the soil-to-soil recovery coefficient of one degree term, static friction factor, rolling friction factor, and soil JKR surface energy and the recovery coefficient of quadratic term soil-to-soil had a highly significant effect on the area difference and the interaction term soil-to-soil static friction factor with soil JKR surface energy had a significant effect on the area difference.

The optimal combination of parameters was determined using the response surface optimization method with the minimum area difference as the target value for the area difference regression model, which comprised the soil-to-soil recovery coefficient of 0.5, the static friction factor of 0.2, the rolling friction factor of 0.05, and the soil JKR surface energy of 5 J/m². At this moment, the relative error of the cross-sectional area of the furrow between the simulation test and the soil tank test was 9.15%, less than the acceptable error of 15% in general statistics and within a reasonable range. The results reveal that the model based on the area difference method can accurately reflect the movement law of soil particles. Moreover, the result of the traction resistance verification test indicated that the relative error between simulation test and soil tank test was 10.4%, which was within a reasonable range. The above results reveal that the mechanical properties of the simulated soil model are basically consistent with the actual soil; it can be applied to the design, selection, use and commissioning of soil contact components in agriculture, and provides technical guidance for the development and design of new soil contact components with low power consumption and ideal operating results.