Modeling and Parameter Selection of the Corn Straw–Soil Composite Model Based on the DEM

Abstract

During corn harvesting operations, machine–straw–soil contact often occurs, but there is a lack of research related to the role of straw–soil contact. Therefore, in this study, a composite contact model of corn straw‒soil particles was established based on the discrete element method (DEM). First, the discrete element Hertz‒Mindlin method with bonding particle contact was used to establish a numerical model of the double-bonded bimodal distribution of corn straw, and bonding particle models of the outer skin‒outer skin, inner pulp‒inner pulp, and outer skin‒inner pulp were developed. The nonhomogeneous and deformable material properties were accurately expressed. The straw compression test combined with simulation calibration was used to determine some of the bonding contact parameters by means of the PB (Plackett–Burman) test, the steepest ascent test, and the BB (Box–Behnken) test. Additionally, Additionally, the Hertz-Mindlin with JKR (Johnson-Kendall-Roberts) + bonding key model was used to establish the numerical model of the soil particles, which was used to describe the irregularity and adhesion properties of the soil particles. The geometric model of the soil particles was established using the multisphere filling method. Finally, a composite contact model of corn straw‒soil particles was established, the contact parameters between straw and soil were calibrated via collision tests, inclined tests and inclined rolling tests, and the established composite contact model was further verified through direct shear tests between straw and soil. A theoretical foundation for the optimal design of equipment linked to maize harvesting is provided by this work.

1. Introduction

During corn harvesting operations, there is often machine–straw–soil contact, so it is necessary to clarify the role of this machine–straw–soil contact. With the development of numerical simulation theory, an increasing number of studies have used the numerical simulation software EDEM (EDEM 2018, DEM Solutions, Edinburgh, UK, 2002) to analyze straw‒machine and machine‒soil interactions [1,2,3,4,5], but few reports on the contact interactions between straw and soil exist. Thus, it is essential to use advanced simulation techniques to improve the understanding of how straw interacts with the soil. Therefore, this research can lead to better machinery designs that enhance efficiency and reduce soil disturbance.

In addition, when modeling straw with a DEM, it is often assumed to be homogeneous and treated as a single medium. For example, Zeng et al. assumed the straw to be a rigid body in the simulation of the straw clearing process, ignored the effect of straw deformation, and simplified it into a string of overlapping spheres of cylinders [6]. Xu et al. simplified and approximated the cucumber straw as a cylinder and built a model of cucumber straw particles using a number of spherical particles in the EDEM units to model cucumber straw particles [7]. Zheng et al. used the DEM to establish a corn straw cutting model, and their straw model was also built from cylindrical filled particles [8]. However, these studies ignored the difference in mechanical properties between the outer skin and inner pulp, which is contrary to the actual situation of straw. Therefore, when a model for the discrete element analysis of corn straw is established, the outer skin and inner pulp should be differentiated to reflect the differences in the physical-mechanical parameters and to improve the modeling accuracy and reliability. Accordingly, Liu Y C et al. established a double-bonded model of corn straw, which was used to represent the mechanical difference between the outer skin and inner pulp, and theoretically calculated the bonding parameters in the bonding model [9]. Liu W et al. established a bonding model of corn straw with different structures and a model of particles with different properties, employing a combination of simulation optimization design techniques and physical testing to validate the discrete unit model; the bonding parameters of the model were thereby calibrated [10]. However, it is essential to conduct an in-depth study on how to increase the filling density and reduce the porosity of the model; this facilitates stronger and more reliable adhesion among particles, thereby meeting the simulation accuracy requirements. Moreover, soil particles are characterized by significant nonhomogeneity and nonlinearity, resulting in a more complex ontological relationship. Their characteristics are significantly influenced by such elements as texture, porosity, organic matter content, and water content. When modeling soil particles with DEMs, the multisphere filling method is generally used to model their geometry, but the geometric model of the soil established in the past is too simple and insufficient to express their complex irregular shape. For example, Xu et al. used the multisphere method and particle size enlargement method to establish a soil particle model based on the number of different filling spheres [11]. In the selection of mechanical models with soil particles, scholars often choose the Hertz–Mindlin model combined with the JKR model to describe the adhesion between different soil particles [12,13,14,15]. The above studies show that the Hertz–Mindlin model combined with the JKR model can describe soil–soil particle adhesion, but its simulation results have some errors with respect to the experimental results that must be corrected further [16]. Therefore, the contact mechanics model used to accurately express soil‒soil adhesion must be analyzed in depth. Finally, we aimed to explore methods for obtaining the contact parameters between straw and soil, which are challenging to acquire experimentally, and to conduct a thorough investigation into how these parameters can be utilized in validating the established straw‒soil composite contact mechanics model.

Therefore, to establish an accurate straw‒soil composite discrete element model, the objectives of this paper are as follows: (1) accurately determine the physical-mechanical parameters of corn straw by adopting the Hertz‒Mindlin model with bonding particle contact to establish a numerical model of the double-bonded bimodal distribution of corn straw, and by using the joint filling of large and small ball particles to improve the filling density and reduce the porosity of the model, thereby increasing the accuracy of the simulation; (2) Determine the complex bonding contact parameters of straw via experiments combined with calibration; (3) Measure the soil properties in the corn planting area, using the Hertz‒Mindlin model combined with the JKR–bonding key model combination to accurately express the adhesion properties of the soil and establish a discrete element model of the soil; (4) Establish the composite straw‒soil contact simulation model, calibrate the straw‒soil contact parameters via collision tests and inclined rolling tests, and verify the composite model using straight shear tests between straw and soil.

2. Materials and Methods

2.1. Discrete Elemental Modeling of Corn Straw

2.1.1. Modeling the Double-Bonded Bimodal Distribution of Corn Straw

To characterize the mechanical differences between the outer skin xylem and the fibrous part of the inner pulp of corn straw, as well as to describe the nonhomogeneity of the straw, two types of bulk particles with different structures and parameters were established via EDEM to fill the outer skin and the inner pulp of the straw, respectively. The geometrical model of the outer skin of the straw was developed using the particle arrangement method [17], which employs subspheres of uniform size with a radius of 0.75 mm. Subsequently, superimposed spheres were utilized to complete the model, which represented the outer skin of the straw. The inner pulp used the particle replacement method of particle factories in EDEM, which was clarified by several preliminary experiments, to build two different sizes of spherical particles with the same parameters in a cylindrical factory. In this process, the large particles (radius of 1 mm) were filled first and the small particles (radius of 0.5 mm) were filled later with the aim of obtaining a more densely packed arrangement structure. Finally, the inner pulp‒inner pulp, outer skin‒outer skin, and outer skin‒inner pulp together were bonded together using the Hertz–Mindlin bonding particle contact method, from which a numerical model of the double-bonded bimodal distribution of corn straw was established, as shown in Figure 1 below.

2.1.2. Selection of Corn Straw Contact Parameters

Corn straw from the experimental field of Jilin Agricultural University that was without obvious breakage or pest symptoms was selected as the test material. After the leaves of the selected straw were removed, only the culm portion was retained for the physical-mechanical test. The moisture content was 9.5 ± 1.5%. The density was determined via the drainage method, and the densities of the inner and outer skin of the straw were 185 kg/m³ and 614 kg/m³, respectively.

The Poisson’s ratio and modulus of elasticity of the straw, inner pulp, and outer skin were measured via uniaxial compression tests in a universal testing apparatus, as seen in Figure 2 below. Each group underwent five repetitions of the single-factor test with a loading rate of 2 mm/min, and the final value was determined by averaging the results. Poisson’s ratio was calculated via Equation (1):

$$\mu =\left|{\displaystyle \frac{{\delta }_1}{{\delta }_2}}\right|={\displaystyle \frac{W_1-W_2}{L_1-L_2}}$$

(1)

where δ₁ is the transverse deformation, mm; δ₂ is the axial deformation, mm; W₁ is the transverse dimension before compression, mm; W₂ is the transverse dimension after compression, mm; L₁ is the axial dimension before compression, mm; and L₂ is the axial dimension after compression, mm.

Finally, the inner pulp Poisson’s ratio of the straw was measured to be 0.26, and Poisson’s ratio for the outer skin was 0.10.

Equation (2) was used to get the modulus of elasticity:

$$E={\displaystyle \frac{FL}{S\Delta L}}$$

(2)

where F is the maximum bearing force, N; L is the initial length, mm; S is the sample cross-sectional area, mm; and ∆L is the difference in length before and after sample compression, mm.

Finally, the moduli of elasticity of the inner pulp and outer skin of the straw were 48.44 MPa and 163.39 MPa, respectively.

Based on the experimental results and some parameters taken from references [10,18], the resulting material intrinsic parameters and contact parameters are shown in

Table 1. Material eigen parameters and contact parameters.

Materials	Parameters	Values	Source
Straw	Density (Kg/m3)	308	Measured
	Poisson’s ratio	0.39	Measured
	Shear modulus (Pa)	4.59 × 107	Measured
Inner pulp	Density (Kg/m3)	185	Measured
	Poisson’s ratio	0.26	Measured
	Shear modulus (Pa)	1.92 × 107	Measured
Outer skin	Density (Kg/m3)	614	Measured
	Poisson’s ratio	0.10	Measured
	Shear modulus (Pa)	7.43 × 107	Measured
Steel	Density (Kg/m3)	7800	Reference [18]
	Poisson’s ratio	0.30	Reference [18]
	Shear modulus (Pa)	7.10 × 1010	Reference [18]
Outer skin-Steel	Collision recovery coefficient	0.608	Reference [10]
	Coefficient of static friction	0.328	Reference [10]
	Coefficient of rolling friction	0.112	Reference [10]
Inner pulp-Steel	Collision recovery coefficient	0.358	Reference [10]
	Coefficient of static friction	0.378	Reference [10]
	Coefficient of rolling friction	0.128	Reference [10]
Outer skin-Outer skin	Collision recovery coefficient	0.481	Reference [10]
	Coefficient of static friction	0.211	Reference [10]
	Coefficient of rolling friction	0.096	Reference [10]
Inner pulp-Inner pulp	Collision recovery coefficient	0.301	Reference [10]
	Coefficient of static friction	0.487	Reference [10]
	Coefficient of rolling friction	0.142	Reference [10]
Outer skin-Inner pulp	Collision recovery coefficient	0.385	Reference [10]
	Coefficient of static friction	0.432	Reference [10]
	Coefficient of rolling friction	0.118	Reference [10]

.

Through the axial compression test of straw combined with the calibration method, the mechanical parameters of the parallel bonding model (PBM) of straw that were difficult to obtain experimentally, such as the normal critical stress, shear critical stress, normal stiffness coefficient, and shear stiffness coefficient, were calibrated, and the optimal parameter combinations were finally obtained.

For the simulation and simulation analysis, the Hertz‒Mindlin model with bonding particle contact was used, and the significance of the simulation parameters was analyzed via the PB test. Then, the levels of the test factors gradually approached the optimal ranges via the steepest ascent test. Finally, the BB test was carried out to derive the optimal combinations of the simulation parameters for the double-bonded bimodal distribution model of corn straw.

First, the literature was reviewed in order to choose the range of bonding parameter values [10,19,20]; next, the range of bonding contact parameters was selected according to the moisture content of the corresponding straw, and the PB test was carried out to screen out the factors that were significant to the test results. The Design Expert 13 software was used, and the discrete element contact parameters of straw were taken as factors, denoted by A~D (the subscripts 1, 2, and 3 represent the contact parameters of the outer skin–outer skin, outer skin–inner pulp, and inner pulp–inner pulp, respectively), in which N~T are virtual parameters. The test values of each factor were used as the level values, and the response value was the compression force.

Table 2. Plackett–Burman test parameter levels.

Variety	Factor	Low Level (−1)	High Level (+1)
A1	Normal stiffness per unit area (N/m3)	4.00 × 109	8.00 × 109
B1	Shear stiffness per unit area (N/m3)	3.00 × 108	9.00 × 108
C1	Critical normal stress (Pa)	1.00 × 108	5.40 × 108
D1	Critical shear stress (Pa)	1.00 × 107	5.00 × 107
A2	Normal stiffness per unit area (N/m3)	2.00 × 108	1.00 × 109
B2	Shear stiffness per unit area (N/m3)	2.00 × 108	6.00 × 108
C2	Critical normal stress (Pa)	1.50 × 107	5.50 × 107
D2	Critical shear stress (Pa)	2.50 × 105	6.50 × 105
A3	Normal stiffness per unit area (N/m3)	2.00 × 107	6.50 × 108
B3	Shear stiffness per unit area (N/m3)	6.00 × 106	6.00 × 107
C3	Critical normal stress (Pa)	4.50 × 105	6.50 × 105
D3	Critical shear stress (Pa)	2.50 × 105	4.50 × 105
N, O, P, Q, R, S, T	Virtual parameters

provides the table of test parameter levels.

The loading speed of the axial simulation test plate was set to 2 mm/min, which was the same running speed as that of the mechanical test. A compression test simulation screenshot is shown in Figure 3 below.

Figure 4 shows the comparison curve graphs of the compression test results of the simulation and the actual experiment. The figure shows that both have the same trend: approximate yielding occurs at the compression displacement near 2 mm. The experimental curve graph reveals that some yielding occurs at a compression displacement of 1.3 mm, while the yield is transitional rather than developmental; this is because parts of the stalk twist during compression, but the stalk may continue to compress when pressure is applied continuously.

The PB test consists of 20 groups of tests, for which the specific arrangements and results are listed in

Table 3. Plackett–Burman test design and results.

NO.	A1	B1	C1	D1	A2	B2	C2	D2	A3	B3	C3	D3	Compression Force (N)
1	+1	+1	−1	−1	+1	+1	+1	+1	−1	+1	−1	+1	157.49
2	−1	+1	+1	−1	−1	+1	+1	+1	+1	−1	+1	−1	111.09
3	+1	−1	+1	+1	−1	−1	+1	+1	+1	+1	−1	+1	166.17
4	+1	+1	−1	+1	+1	−1	−1	+1	+1	+1	+1	−1	201.18
5	−1	+1	+1	−1	+1	+1	−1	−1	+1	+1	+1	+1	152.05
6	−1	−1	+1	+1	−1	+1	+1	−1	−1	+1	+1	+1	90.63
7	−1	−1	−1	+1	+1	−1	+1	+1	−1	−1	+1	+1	53.57
8	−1	−1	−1	−1	+1	+1	−1	+1	+1	−1	−1	+1	136.71
9	+1	−1	−1	−1	−1	+1	+1	−1	+1	+1	−1	−1	168.89
10	−1	+1	−1	−1	−1	−1	+1	+1	−1	+1	+1	−1	109.39
11	+1	−1	+1	−1	−1	−1	−1	+1	+1	−1	+1	+1	128.32
12	−1	+1	−1	+1	−1	−1	−1	−1	+1	+1	−1	+1	159.80
13	+1	−1	+1	−1	+1	−1	−1	−1	−1	+1	+1	−1	102.93
14	+1	+1	−1	+1	−1	+1	−1	−1	−1	−1	+1	+1	83.32
15	+1	+1	+1	−1	+1	−1	+1	−1	−1	−1	−1	+1	90.12
16	+1	+1	+1	+1	−1	+1	−1	+1	−1	−1	−1	−1	98.04
17	−1	+1	+1	+1	+1	−1	+1	−1	+1	−1	−1	−1	127.08
18	−1	−1	+1	+1	+1	+1	−1	+1	−1	+1	−1	−1	119.45
19	+1	−1	−1	+1	+1	+1	+1	−1	+1	−1	+1	−1	125.26
20	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	87.34

below.

Table 4. Analysis of the Plackett–Burman test factors’ significance.

Factor	Sum of Squares	F Value	p Value	Significance
Model	22,147.82	10.10	0.0026	significance
A1	1195.21	6.54	0.0377	3
B1	848.77	4.65	0.0681	5
C1	296.99	1.63	0.2430	9
D1	0.0014	7.91 × 10−6	0.9978	12
A2	343.21	1.88	0.2128	8
B2	0.4410	0.0024	0.9622	11
C2	400.07	2.19	0.1825	7
D2	649.69	3.56	0.1013	6
A3	10,777.33	58.99	0.0001	1
B3	6739.22	36.89	0.0005	2
C3	889.11	4.87	0.0632	4
D3	7.78	0.0426	0.8424	10

shows the significance analysis of the factors tested in the PB test, and the results of the test show that the model has a p value of 0.0026 < 0.01, which is highly significant, making the model reliable; the coefficient of determination R² = 0.9454 indicates that the model is well correlated, has a good fit to the experimental data, and is applicable to 94.54% of the experimental data; the adjusted coefficient of determination, adjusted R²_adj = 0.8518, indicates that the model explains 85.18% of the response variable’s variability after accounting for the quantity of independent variables in the model; p < 0.01 for A₃ and B₃, indicating that normal stiffness per unit area of inner pulp–inner pulp and shear stiffness per unit area of inner pulp–inner pulp have a highly important impact on the critical load; 0.01 < p < 0.05 for A₁, indicating that normal stiffness per unit area of outer skin–outer skin has a significant effect on the critical load; p > 0.05 for the other factors, suggesting that there is no discernible impact on the critical load.

Combined with Figure 5, the Pareto chart for the PB test shows that factors A₁, A₃, and B₃ all exceeded the t = 2.36462 value and all had a favorable impact on the PB test response value, with A₃ being the most significant factor, indicating the direction of ascent in the steepest ascent test path.

In summary, the steepest ascent test was conducted for the three most significant factors—A₁, A₃, and B₃.

Table 5. Arrangement and outcomes of the steepest ascent method’s path.

No.	1	2	3	4	5
A1	4.00 × 109	5.00 × 109	6.00 × 109	7.00 × 109	8.00 × 109
A3	2.00 × 107	1.775 × 108	3.35 × 108	4.925 × 108	6.50 × 108
B3	6.00 × 106	1.95 × 107	3.30 × 107	4.65 × 107	6.00 × 107
Compression force(N)	70.13	82.96	95.00	105.11	111.43
Relative Error	39.4%	28.3%	17.8%	9.1%	3.6%

displays the precise path configurations and outcomes. An examination of the results of the steepest ascent test, as shown in Table 5, revealed that the critical load at the fifth set of test levels was 111.43 N, with the smallest relative error of 3.6% compared with the critical load of 115.64 N in the physical test. Here is observed the deviation between the measured value and the real value of the critical load on the corn straw. The relative error was calculated via Equation (3):

$$\delta ={\displaystyle \frac{\left|F_P-F_S\right|}{F_P}}\times 100\%$$

(3)

where F_P is the measured value of compressive force, N, and F_S is the simulation value of the compression force, N.

However, for the fifth group of tests as an endpoint and the critical load being smaller than the critical load of the actual experiment, the choice of high level should be adjusted. At this point, the sixth set of tests was added, and when A₁ was 9.00 × 10⁹ N/m³, A₃ was 8.075 × 10⁸ N/m³, B₃ was 7.35 × 10⁷ N/m³, and the critical load was 117.08 N, with a relative error of 1.2%. Therefore, Test 5 was selected as the center point in the follow-up test, and the factor levels of Test 4 and Test 6 were used as the low and high levels of the BB test, respectively, for the response surface test with three factors and three levels analysis.

Table 6. Design and outcomes of the Box–Behnken test.

No.	Factor Level Value			Compression Force (N)
	A1	A3	B3
1	−1	−1	0	105.93
2	+1	−1	0	110.94
3	−1	+1	0	115.6
4	+1	+1	0	107.07
5	−1	0	−1	111.42
6	+1	0	−1	105.27
7	−1	0	+1	109.36
8	+1	0	+1	112.61
9	0	−1	−1	104.28
10	0	+1	−1	109.16
11	0	−1	+1	108.29
12	0	+1	+1	111.5
13	0	0	0	115.64
14	0	0	0	116.51
15	0	0	0	116.42
16	0	0	0	115.62
17	0	0	0	113.74

shows the BB test design and results. In this case, B₁ was 6.00 × 10⁸ N/m³, C₁ was 3.20 × 10⁸ Pa, and D₁ was 3.00 × 10⁷ Pa; A₂ was 6.00 × 10⁸ N/m³, B₂ was 4.00 × 10⁸ N/m³, C₂ was 3.00 × 10⁷ Pa, and D₂ was 4.50 × 10⁵ Pa; C₃ was 5.5 × 10⁵ Pa, and D₃ was 3.50 × 10⁵ Pa.

By fitting a quadratic multiple regression to the experimental results, the regression equation with the critical load as the response value was obtained as follows:

Y = 115.59 − 0.8025A₁ + 1.74A₃ + 1.45A₃ − 3.39A₁A₃ + 2.35A₁B₃ − 0.4175A₃B₃ − 2.17A₁² − 3.53A₃² − 3.75B₃². The coefficient of determination is R² = 0.9782, the adjusted coefficient of determination adjusted R²adj = 0.9502, and the coefficient of variation is C.V. = 0.8196%. The regression equation indicated that the model fit was high enough to explain 95.02% of the variability in the response values, so ANOVA was performed on the BB test, the results of which are shown in

Table 7. Significance analysis using the Box–Behnken test.

Source	Sum of Squares	df	Mean Square	F-Value	p-Value
Model	260.73	9	28.97	34.92	<0.0001	significant
A1	5.15	1	5.15	6.21	0.0415
A3	24.12	1	24.12	29.07	0.0010
B3	16.91	1	16.91	20.38	0.0028
A1A3	45.83	1	45.83	55.24	0.0001
A1B3	22.09	1	22.09	26.62	0.0013
A3B3	0.6972	1	0.6972	0.8403	0.3898
A12	19.86	1	19.86	23.94	0.0018
A32	52.44	1	52.44	63.21	<0.0001
B32	59.19	1	59.19	71.34	<0.0001
Residual	5.81	7	0.8297
Lack of Fit	0.8467	3	0.2822	0.2275	0.8730	not significant
Pure error	4.96	4	1.24
Cor Total	266.53	16

. The model has p < 0.01, indicating that the model is highly significant and can reflect the connection between response values and factors—whereas the misfit term is not significant (p > 0.05)—and that there is no significant misfit problem. The primary factors A₃ and B₃ have highly significant effects on the critical loads, and the effect of A₁ is significant. The effect of A₃B₃ is negligible, but the interaction terms A₁A₃ and A₁B₃ are considerable, and the effects of secondary terms A₁², A₃², and B₃² are highly significant. The effects of the factors on the response values did not have a simple linear relationship in the order of A₁ < B₃ < A_3.

The regression equation was used to predict the compression force for each bonding parameter condition, and the relationship between the physical and simulation tests of the model was plotted (Figure 6), as were the contour maps (Figure 7). Based on the correlation between the expected and experimental values (Figure 6), the distribution of the data points was close to the 45° fitting line and was more concentrated, which indicated that the model fit better and was more accurate in terms of prediction, and thus that it could be used for predicting and analyzing the compression force of corn straw. The contour maps show the trends and interactions of the effects of the factors, as well as the extent of the experiment conducted. The response variable that reached the optimal value appears in the red region. Combined with the BB test analysis of each optimal parameter value, A₁ was 8.20 × 10⁹ N/m³, A₃ was 6.81 × 10⁸ N/m³, and B₃ was 6.34 × 10⁷ N/m³.

In addition, the contact radius of the Hertz‒Mindlin bonding contact mechanics model is generally 1.2‒2 times the physical radius of the particles [19]. To avoid bonding between non-neighboring particles or the instability of bonds formed by neighboring particles, many preliminary experiments concluded that the particles’ interaction radius was set to 1.2 times their own radius, which not only had a better simulation effect but also enabled the corn straw model to have certain flexible characteristics. Therefore, the outer skin particles had a contact radius of 0.9 mm, the inner pulp big particles had a contact radius of 1.2 mm, and the inner pulp small particles had a contact radius of 0.6 mm. The bonding radius and the contact radius were taken to be the same value, where the physical radius of dissimilar particles was taken as 1.2 times the center value, and the outer skin particle–inner pulp large particle bonding radius was 1.05 mm. The outer skin particle–inner pulp small particle bonding radius was 0.75 mm, and the inner pulp large particle–inner pulp small particle bonding radius was 0.9 mm.

2.2. Discrete Elemental Modeling of Soils

2.2.1. Modeling Soil Particles Hertz–Mindlin with JKR + Bonding

To describe the diversity of soil particle shapes, the soil particles were modeled in the form of multisphere assemblies. Nonregular particles were modeled using the discrete element software EDEM to create soil particles similar to spheres, cylinders, pie-shaped particles, triangles, and quadrangular prisms with radial dimensions of 2.2–3.0 mm, as shown in Figure 8 below.

For sandy loam soils with water contents of 20% and above, owing to the existence of adhesion between soil particles, it was necessary to consider interparticle adhesion, and the normal and shear forces between the contacting particles were determined using the Hertz-Mindlin with JKR+ bonding model [21,22].

2.2.2. Selection of Soil Particle Contact Parameters

In accordance with the preliminary work of the research group, a combination of simulations and experiments was used to calibrate the contact parameters between the soils via the steepest ascent test and BBD (Box‒Behnken Design) test combined with the direct shear test. The soil parameters were obtained from the experiment and BBD test analysis, with a density of 1950 kg/m³, shear modulus of 2.73 × 10⁶ Pa, Poisson’s ratio of 0.2, coefficient of restitution of 0.3, coefficient of static friction of 0.5, coefficient of rolling friction of 0.03, surface energy of 4.436 J/m², normal stiffness per unit area of 2.86 × 10⁶ N/m³, shear stiffness per unit area of 5.54 × 10⁵ N/m³, critical normal stress of 1833 Pa, and critical shear stress of 3332 Pa [23].

2.3. Modeling and Validation of a Corn Straw–Soil Composite Contact Model

To obtain more accurate straw‒soil contact parameters, the collision recovery coefficient, static friction coefficient, and rolling friction coefficient of straw‒soil particles were calibrated using collision tests, inclined tests, and inclined rolling tests, respectively, in the form of a combination of experimental and simulation tests via a single-factor test.

The collision recovery coefficient of straw–soil was calibrated using the collision test. A high-speed camera was used to record the highest height of the straw after free fall from a height of 200 mm to the soil bin. Five repeated tests were conducted, and the average value of 0.26 was taken as the value of the straw‒soil collision recovery coefficient, as shown in Figure 9. In the EDEM simulation test, the static friction factor and rolling friction factor of straw–soil were considered to be 0. After the preliminary simulation tests, the collision recovery coefficient of straw–soil was set to 0.15, 0.2, 0.25, 0.3, and 0.35 to carry out five sets of simulation tests; the simulation process is shown in Figure 9 below. The simulation results are shown in Figure 12a below, indicating that as the collision recovery coefficient rises, the results of the simulation varied back and forth around the experimental mean value within the error band; moreover, when the collision recovery coefficient of straw–soil was set to 0.2, the simulation result was closest to the experimental value, so the collision recovery coefficient of straw and soil was 0.2.

The static friction coefficient of straw–soil was calibrated via the inclined test. The straw was put on a horizontal groove of soil and, after rotating the groove steadily and slowly within the fixed axis, the rotation was halted when the straw started to slip. The angle of inclination of the rotating plate was recorded as shown in Figure 10 below, and the coefficient of static friction between the straw and the soil was 0.21 after five repetitions of the experiment to obtain the average value. In the EDEM simulation test, the rolling friction coefficient value of straw–soil was taken as zero, whereas the collision recovery coefficient of straw–soil was adopted as the calibrated value of 0.2. After the preliminary simulation tests, the static friction coefficient of straw–soil was set to 0.15, 0.20, 0.25, 0.30, and 0.35 for the five sets of simulation tests, the simulation screenshot of which is shown in Figure 10 below. The simulation results are shown in Figure 12b below, indicating that with the exception of the straw–soil static friction coefficient, which was 0.20, the simulation results varied as the static friction coefficient increased. The remaining values fell within the error band, and the simulation result was closest to the experimental value when the straw–soil static friction coefficient was 0.25. Therefore, it was determined that the straw–soil static friction coefficient was 0.25.

The rolling friction coefficient of straw–soil was calibrated via the inclined rolling test, which is shown in Figure 11. The straw was placed on an inclined soil groove with an inclination angle of 30° and a fixed height of 7.0 cm, and it rolled down the sloping earth groove from a starting speed of zero. The straw rolled down the horizontal surface until it reached a complete standstill, and the horizontal rolling distance obtained was measured. Five repetitions of the experiment were carried out to obtain an average value of 8.38 cm. In the EDEM simulation test, the collision recovery factor and static friction factor of straw–soil were used, which were calibrated at 0.2 and 0.4, respectively. After the preliminary simulation tests, the rolling friction coefficient between straw and soil was set to 0.08, 0.1, 0.12, 0.14, and 0.16, and five sets of simulation tests were carried out. A simulation screenshot is shown in Figure 11. The simulation results are shown in Figure 12c, indicating that the simulation results fluctuated back and forth with increasing rolling coefficient; the simulation result was closest to the experimental value when the rolling friction coefficient of soil and straw was 0.1.

2.3.1. Corn Straw‒Soil Composite Direct Shear Test

The straw‒soil composite direct shear test was carried out using a ZJ-type strain-controlled direct shear apparatus, as shown in Figure 13. The shear box was evenly filled with the soil–straw mixture and left with a smooth surface. Normal compressive loads of 50 kPa, 100 kPa, 150 kPa, and 200 kPa were applied to each sample. For every vertical pressure, the experimental data were recorded to obtain the shear strength of the soil–straw mixtures as a function of vertical loading under different vertical stresses.

2.3.2. Simulation of the Corn Straw‒Soil Composite Direct Shear Test

When simulating the straw‒soil composite direct shear test process, the bonding key was included, and the soil particles were simulated using the JKR model; moreover, a numerical model of the double-bonded bimodal distribution of corn straw was developed by combining the discrete element Hertz‒Mindlin model with the bonding particle contact method. To guarantee the stability and convergence of the simulation’s numerical computations, the simulation’s time step was set to 5 × 10⁻⁷ s, and the simulation duration was 9.5 s.

The loading speed of the shear box under the straw‒soil composite direct shear simulation test was set to 0.8 mm/min to maintain the same speed as the running speed of the mechanical test; the screenshot of the shear test simulation is shown in Figure 14 below.

3. Results

A comparison between the simulation and experimental results is presented in Figure 15. The figure shows that the shear strength of the straw–soil composites increased with increasing normal loads in both the simulations and tests. Moreover, the simulation results showed the same trend as the experimental data, and the deformation of the corn straw remained largely unchanged before and after the testing process, thereby verifying the feasibility of the model established in this study.

Furthermore, to better observe the deformation of the numerical model of the double-bonded bimodal distribution of straw during shearing, all the soil particles were hidden and the deformation maps of the straw in the +x and +z directions were intercepted, as shown in Figure 16. The figure shows that with the application of a normal load, the interaction between the straw and the soil, as well as the relative movement of the upper and lower shear boxes, significantly changed the deformation and displacement of the straw.

To better observe the deformation of straw under different normal loads in the numerical model of the double-bonded bimodal distribution of straw during the shearing process, all the soil particles were hidden, and the deformation map of straw in the +z direction was intercepted, as shown in Figure 17 below. The figure clearly shows that the deformation of the straw varies significantly with the application of a normal load, as well as the interaction between the straw and the soil, the relative movement of the upper and lower shear boxes, and the fact that the deformation of the straw increases with increasing normal load. Moreover, the number of internal bonding key ruptures increased. The discrete element model of straw without an applied load produced a total of 550 intact bonding keys, which gradually ruptured as the load was applied and increased; the number of bonding keys gradually decreased with increasing normal load. When the load reached 50 kPa, the number of intact bonds was 503. When loaded to 100 kPa, the number of intact bonds was 465. When the load reached 150 kPa, the number of intact bonds was 454, and when the load reached 200 kPa, the number of intact bonds was only 390.

4. Discussion

In this study, corn straw–soil was taken as the research object, and a straw discrete element model, soil discrete element model and straw–soil composite discrete element model were sequentially established via the DEM.

4.1. Particle Modeling

The corn straw was treated in layers. The inner pulp of the straw was filled with large and small ball particles of two different sizes, and the outer skin was modeled using the particle arrangement method. The non-homogeneous and deformable material properties were accurately expressed, thereby avoiding the limitations of previous studies by Zeng and Xu et al. [6,7,8,9,10] and clearly indicating the difference in mechanical properties between the outer skin and the inner pulp. The filling density can be improved and the reduced porosity of the model can increase the accuracy of the simulation.

4.2. Contact Model Selection

The corn straw model adopted the Hertz–Mindlin combined with the bonding model. This model enables adjacent particles to have parallel bonding at the contact points, facilitating the particles in acquiring a more stable bonding force and enhancing the precision of the model. When exposed to external forces, the bonding keys will fracture, and the related mechanical characteristics, such as the strength of the model, can be reflected through the rupture of the bonding keys. The Hertz–Mindlin combined with the JKR– + bonding model combination was used to describe the adhesion characteristics of soils. The contact models were sufficient to describe the complicated dynamic behavior between straw and soil. These conclusions are consistent with Liu Y.C and Min Liu et al. [9,13].

4.3. Model Accuracy Verification

There are few reports on straw–soil contact interactions. It is difficult to find a suitable validation method for the straw–soil composite model. The previous verification methods include the pilling test, the bulk density test, the rotated drum test, and the “self-screening” test [24,25,26,27,28,29,30,31]. In this paper, the direct shear test for straw–soil composites was first adopted for model verification. According to the simulation results, this approach can improve the understanding of the interaction between straw and soils.

5. Conclusions

In this work, corn straw–soil was taken as the research object, and the straw discrete element model, the soil discrete element model, and the straw–soil composite discrete element model were sequentially established via the DEM. Moreover, the feasibility and effectiveness of the straw‒soil composite modeling method proposed in this paper were verified through simulations and experimental comparisons of straw−soil composite direct shear tests, and the conclusions are as follows:

• The particle arrangement method was used; when modeling straw, the outer skin and inner pulp were layered. A numerical model of the double-bonded bimodal distribution of corn straw was established via the Hertz‒Mindlin bonding contact model, and the outer skin and inner pulp of the straw were filled with two different sizes of large-ball particles and small-ball particles. Soil particles were modeled in the form of multisphere combinations, and soil particles were modeled via the Hertz‒Mindlin with JKR+ Bonding model.
• The axial compression test of straw combined with the calibration method was used to determine the mechanical parameters of the parallel bond model (PBM) of straw, which are difficult to obtain through experimentation. The factors that had a significant effect on the critical load of axial compression were screened via the PB (Plackett–Burman) test. Afterwards, the steepest ascent test was conducted to bring the test factor levels to the optimal range, and finally, the BB (Box‒Behnken) test was conducted to derive the optimal simulation parameter combinations: outer skin–outer skin normal stiffness per unit area of 8.20 × 10⁹ N/m³, inner pulp–inner pulp normal stiffness per unit area of 6.81 × 10⁸ N/m³, and inner pulp–inner pulp shear stiffness per unit area of 6.34 × 10⁷ N/m³.
• A straw‒soil composite discrete element model was established. The exposure parameters of the straw‒soil composite model were calibrated via three single-factor tests, including collision tests, inclined tests and inclined rolling tests; the collision recovery coefficients, static friction coefficients, and rolling friction coefficients of straw–soil were 0.2, 0.25, and 0.1, respectively.
• To validate the model, a straw–soil composite direct shear test was employed. The results showed that the simulation results had the same trend as the experimental results, and that the values were close to each other. The reliability and reasonableness of the model parameter selection and modeling method proposed were confirmed, and a reasonable straw‒soil composite discrete element model and simulation parameters were provided for the simulation and optimization design of maize harvesting operation-related machinery.