Calibration and Verification of Discrete Element Parameters of Surface Soil in Camellia Oleifera Forest

Abstract

To analyze the interaction between the surface soil and the soil-contacting component (65 Mn) in the camellia oleifera forest planting area in Changsha City, Hunan, China, in this study, we conducted discrete element calibration using physical and simulation tests. The chosen contact model was Hertz–Mindlin with JKR cohesion, with the soil repose angle as the response variable. The repose angle of the soil was determined to be 36.03° based on the physical tests. The significant influencing factors of the repose angle determined based on the Plackett–Burman test were the soil–soil recovery coefficient, soil–soil rolling friction coefficient, soil-65 Mn static friction coefficient, and surface energy of soil for the JKR model. A regression model for the repose angle was developed using the Box–Behnken response surface optimization method to identify the best parameter combination. The optimal parameter combination for the JKR model was determined as follows: surface energy of soil: 0.400, soil–soil rolling friction coefficient: 0.040, soil-65 Mn static friction coefficient: 0.404, and soil–soil recovery coefficient: 0.522. The calibrated discrete element parameters were validated through experiments on the repose angle and steel rod insertion. The results indicated that the relative errors obtained from the two verification methods were 2.44% and 1.71%, respectively. This research offers fundamental insights for understanding the interaction between soil and soil-contacting components and optimizing their design.

1. Introduction

Camellia oleifera, one of the four major edible raw oil materials worldwide, plays a significant role in China’s agricultural economy. Due to its well-developed root system, drought tolerance, and ability to thrive in less fertile areas, Camellia oleifera is widely distributed across various regions in China. Hunan, in particular, is considered a suitable region for cultivating Camellia oleifera forests. However, the mountainous and hilly terrain of Camellia oleifera plantations makes mechanized operations challenging, resulting in low levels of mechanization and slow industry growth [1,2,3].

The discrete element is a numerical simulation method that can treat the entire medium as a collection of several particle units. It is widely used in the fluidity of scattered materials, solid crushing, and machine–soil interactions [4,5]. Due to the complex soil characteristics, the finite element soil model is not very accurate, and can only simulate soil damage behavior, but cannot simulate the soil movement process. The discrete element method can solve the contact between particles and boundaries and uses adhesive particles to simulate the generation of soil aggregates, greatly improving the accuracy of soil models [6,7].

Rotary tillage plays a crucial role in the management of land operations. It transforms the soil from a compacted state into granular particles with a specific particle size. By studying the interaction between the rotary blade and the soil, valuable insights can be gained to inform blade design. In a study conducted by Cheng et al. [8], the rotary tillage process of a rotary tiller blade was simulated using the discrete element method. The findings indicated that the soil adhesion force primarily accumulates on the oblique long part of the rotary tillage blade. Moreover, as the rotary tillage section moves through the soil, it collects more soil particles. To study the interaction mechanism between vertical rotary tillers and soil, Shike Zhai et al. [9] used EDEM to conduct a virtual simulation of different operating parameters and structural parameters of vertical rotary tillers and obtained the optimal operating parameters and structural parameters. Xiongye Zhang et al. [10] utilized a simulation based on the discrete element method (DEM) to develop the contact model, soil particle model, and soil-rotary tiller roller interaction model. They simulated the dynamic process of the rotary tiller roller cutting soil and obtained information on the soil deformation area, cutting process energy, cutting resistance, and soil particle movement.

Establishing a discrete element model of soil and calibrating contact parameters can enhance the accuracy of numerical simulations in soil and tool preparation. It also provides a foundation and basic parameters for studying soil–tool interaction mechanisms [11,12,13]. Before conducting a DEM simulation analysis, it is crucial to accurately construct the discrete element model of the soil. Yang et al. [14] used the repose angle as a response, calibrated red clay in hilly areas based on physical experiments and simulation tests, and obtained the contact parameters. Du et al. [15] used the Hertz–Mindlin and JKR contact models in EDEM to calibrate the contact parameters involved in the interaction between tilled loose soil and Q235 steel in response to the angle of repose and obtained their optimal parameter combination. The calibration results were verified through field experiments. Zhong et al. [16] used the same method to obtain the discrete element model parameters of soil in rice fields.

Under different scenarios and different operating methods, the soil conditions are different, and the obtained discrete element parameters are also different. Previous research on soil calibration often focuses on using the repose angle as a response variable, with the calibration results validated using the repose angle obtained through physical testing. However, in economic forests, the surface soil tends to be in a fragmented state. When calibrating this type of soil, it is suitable to consider the repose angle as the response variable. Subsequent measurements of the working resistance of soil-contacting components generally align well with the simulation results. Beneath the surface soil lies the root–soil complex, which is no longer in a discrete state due to the significant influence of tree roots. The parameters calibrated using the repose angle as a response show significant errors. In this scenario, it is more suitable to consider the shear stress of the soil or the working resistance of the soil contact parts as a response. To enhance the rigor of the research, this study focused on the surface soil of the Camellia oleifera forest and utilized equal-diameter spherical particles as the soil particle model in order to enhance the simulation’s accuracy.

During the cultivation of Camellia oleifera trees, the interaction between the soil-contacting parts and the soil is involved. Therefore, to address the limited research and lack of references on soil and contact material parameters in the Camellia oleifera forest scene in Changsha, Hunan, this study focused on using local soil as the research object. The Hertz–Mindlin approach with the JKR contact model and the soil repose angle test were employed to calibrate the soil’s discrete element simulation parameters. The calibrated parameters could be utilized for the discrete element simulation between the soil and soil-contacting parts, thereby providing essential data for analyzing the interaction between the soil-contacting parts and soil in Camellia oleifera forest cultivation equipment. Furthermore, it enabled the exploration of the action mechanism of the soil-contacting components and optimization design.

2. Materials and Methods

2.1. Materials

The studied soil was collected on 8 December 2023 from the Camellia oleifera forest experimental base in the Wangcheng District, Changsha City, Hunan Province, China.

In the study, the soil samples were collected from the soil layer of the Camellia oleifera forest at depths ranging from 0 to 50 mm. The sampling method employed was the five-point method, with each sampling point having a minimum sampling mass of 200 g.

The equipment used for the subsequent measurements of the density, friction angle, repose angle, etc., included 200 mm × 400 mm 65 Mn plates, graduated cylinders, balances, cameras, funnels, etc. The software used for these measurements included EDEM and Origin, among others.

2.2. Methods

The Hertz–Mindlin model with JKR was selected as the soil contact model for this study [17]. The parameters, such as the soil’s Poisson’s ratio, solid density, shear modulus, and soil-65 Mn static friction coefficient, were determined through experimental measurements and a literature review. The soil moisture content was obtained using the drying method, the soil-65 Mn static friction coefficient was determined through slope tests, the soil particle size distribution data were obtained through soil sieve screening tests, and the soil density was obtained using an immersion method.

The simulation test was conducted based on EDEM software. The index of the test results was the ratio of the difference between the measured repose angle and the simulated repose angle to the measured repose angle, which was calculated as shown in Equation (1). Taking the above parameters as influencing factors, a second-order regression model of the repose angle and significance parameters was obtained through the Plackett–Burman test, steepest climb test, and Box–Behnken test design, and finally, the best combination of the significance parameters was determined [18,19].

$$y=\frac{\alpha -{\alpha }_1}{{\alpha }_1}\times 100\%$$

(1)

where α₁ is the measured repose angle and α is the discrete element simulation repose angle.

The soil discrete element parameter calibration process is shown in Figure 1 [20,21].

3. Acquisition of Soil Characteristic Parameters and Contact Parameters

3.1. Acquisition of the Soil Characteristic Parameters

The soil characteristic parameters mainly included the moisture content, density, shear modulus, and Poisson’s ratio.

3.1.1. Determination of the Soil Moisture Content

During the winter season, the absence of precipitation persisted for several days, resulting in the soil being dry. The moisture content was determined using the drying method, which revealed that the moisture content of the soil was 1.48%.

3.1.2. Soil Particle Size Distribution Determination

The soil samples were naturally air-dried and weighed using an electronic balance. Subsequently, the samples were sieved using a standard sieve with apertures of 5, 2, 1, 0.5, 0.1, and 0.025 mm [22]. The sieving process was carried out for 10 min (Figure 2).

After sieving, the masses corresponding to different particle sizes were measured using an electronic scale. The mass fractions corresponding to the particle sizes of the soil can be found in

Table 1. Size classification and mass fraction of soil.

Particle size range/mm	0~0.025	0.025~0.1	0.1~0.5	0.5~1	1~2	2~5	>5
Mass fraction/%	0.75	2.91	28.06	18.81	20.62	20.34	8.50

and Figure 3. The mass fractions of the soil particles in different size ranges were as follows: 0.75% (0~0.025 mm), 2.91% (0.025~0.1 mm), 28.06% (0.1~0.5 mm), 18.81% (0.5~1 mm), 20.62% (1~2 mm), 20.34% (2~5 mm), and 8.50% (>5 mm).

3.1.3. Acquisition of Other Soil Characteristic Parameters

Based on previous studies [14,23], the Poisson ratio of the soil was assigned a value of 0.25~0.45, while the shear modulus of the soil was set at 1.15 × 10⁷ Pa.

3.2. Determination of the Soil Static Friction Coefficient

During actual contact between the soil and other materials, static friction plays a major role. The soil static friction coefficient is determined using the inclined plane sliding test [24]. The relevant principles and formulas of its measurement are shown in Figure 4 and Formula (2). The test device comprises two main components: a self-made 200 mm × 400 mm plate made of 65 Mn material and an inclinometer.

To measure the friction angle of soil-65 Mn, we placed the soil sample on the 65 Mn plate and gradually increased the tilt angle of the plate. Once the soil particles on the slope started to slide off, we stopped increasing the slope angle and recorded the tilt angle. This recorded angle value represented the friction angle of soil-65 Mn (Figure 5). We repeated this measurement five times and calculated the average value as the friction angle of soil-65 Mn. The test results are presented in

Table 2. Measurement results of the soil-65 Mn steel static friction angle.

Soil Sample	Soil-65 Mn Static Friction Angle/°
	Test 1	Test 2	Test 3	Test 4	Test 5	Average Value	Overall Average
1	21.40	22.31	22.14	23.22	20.87	22.008	21.9440
2	20.62	21.79	23.41	22.93	20.65	21.88
3	22.82	20.39	21.47	22.83	22.21	21.944

. By using Formula (2), the static friction coefficient of soil-65 Mn was determined as 0.40.

$$\begin{array}{c}{F}_1=G\sin \theta \\ F_2=G\cos \theta \\ F_1-f=0\\ F_1-N=0\\ F_1=\mu F_2\end{array}\}\Rightarrow \mu =\tan \theta$$

(2)

3.3. Determination of the Soil Repose Angle

The repose angle refers to the angle at which a significant amount of material collects on a horizontal surface to form a stable cone-shaped pile. To measure the soil repose angle, we employed the funnel method and constructed a test platform, depicted in Figure 6 [25,26].

In the experiment, the soil was placed into a funnel and allowed to flow from the bottom of the funnel onto the bottom plate, forming a pile shape. A photograph was taken to capture the soil pile, and the photo was subsequently processed to keep only the part to be calculated for the repose angle. After obtaining the soil accumulation photos, image processing technology was used to obtain the repose angle. The general process was as follows:

(1)

We divided the soil pile image into left and right sides symmetrically and extracted the unilateral image;

(2)

We binarized the unilateral image of the soil pile;

(3)

We extracted the image boundary contours;

(4)

The boundary contour line was linearly fitted using the least squares method to obtain the fitting straight line and equation. The slope (k) of the equation was obtained. The repose angle was calculated using Formula (3).

$$\theta =\arctan \left|k\right|$$

(3)

The specific image processing process is shown in Figure 7.

Through the above method, the measurement of the soil repose angle was repeated five times, and the obtained repose angle results are shown in Figure 8.

After the above measurements, it was concluded that the average repose angle of the soil was 36.03°.

4. Soil Discrete Element Parameter Calibration Process and Results

4.1. Construction of the Soil Particle Model and Device Model

The analysis of the soil particle size distribution revealed that the soil particles were predominantly small and exhibited variation in size. Considering the negligible roundness of most soil particles, they could be approximated as standard spherical shapes. To expedite the creation of soil particle models and streamline the subsequent simulation tests, spherical particles with a radius of 1.5 mm were utilized as the soil particle model [27,28,29].

The repose angle detection device model was constructed through SolidWorks software (Figure 9).

4.2. Selection of the Contact Model

The Hertz–Mindlin (no slip) model was chosen as the contact model for the soil-65 Mn interaction in the discrete element simulation, and the Hertz–Mindlin model with JKR cohesion was selected as the contact model for the soil–soil interactions.

The Hertz–Mindlin with JKR cohesion contact model is a cohesive particle contact model that is based on Hertz’s theory. It is specifically designed for simulating materials where particles exhibit bonding and agglomeration due to static electricity, moisture, and other factors. Figure 10 illustrates the schematic diagram of this model. The model accurately represents the viscoelastic characteristics between the particles through normal elastic force, which is determined using the amount of normal overlap and surface energy. The calculation method for normal elastic force is provided in Formulas (4) and (5).

$$F_{JKR}=-4\sqrt{\pi \gamma R^{\ast }{E}^{\ast }}{\delta }_{\tau }^{\frac{3}{2}}+\frac{4E^{\ast }}{3R^{\ast }}{\alpha }^3$$

(4)

$${\delta }_n=\frac{{\delta }_{\tau }^2}{R^{\ast }}-\sqrt{\frac{4\pi \gamma {\delta }_{\tau }}{E^{\ast }}}$$

(5)

In the formula, F_JKR represents the normal elastic force of JKR in Newtons (N), δ_n represents the normal overlap between the two contacting particles in meters (m), δ_τ represents the tangential overlap between the two contacting particles in meters (m), γ represents the surface energy (J/m²), E* represents the equivalent elastic modulus in Pascals (Pa), and R* represents the equivalent radius of the particle in meters (m).

The definitions of the equivalent elastic modulus E* and the equivalent contact radius R* are shown in Formulas (6) and (7).

$$\frac{1}{E^{\ast }}=\frac{1-{\mu }_i^2}{E_i}+\frac{1-{\mu }_j^2}{E_j}$$

(6)

$$\frac{1}{R^{\ast }}=\frac{1}{R_i}+\frac{1}{R_j}$$

(7)

In the formula, E_i and E_j represent the elastic modulus of particle i and particle j in contact (Pa); μ_i and μ_j represent the Poisson ratios of particle i and particle j; and R_i and R_j represent the radii of particle i and particle j (m).

When the surface energy γ = 0, the JKR normal elastic force becomes the Hertz–Mindlin normal force (Formula (8)).

$$F_{JKR}=F_{Hertz}=\frac{4}{3}{E}^{\ast }\sqrt{R^{\ast }}{\delta }_n^{\frac{3}{2}}$$

(8)

When the particles are not in direct contact, the Hertz–Mindlin with JKR contact model can also account for attractive cohesion. Equations (9) and (10) demonstrate the methods to calculate the normal maximum gap and tangential maximum gap when there is non-zero cohesion between the particles.

$${\delta }_c=\frac{{\alpha }_c^2}{R}-\sqrt{\frac{4\pi \gamma {\alpha }_c}{E^{\ast }}}$$

(9)

$${\alpha }_c={\left[\frac{9\pi \gamma R^{\ast 2}}{2E^{\ast }}-\left(\frac{3}{4}-\frac{1}{\sqrt{2}}\right)\right]}^{\frac{1}{3}}$$

(10)

In the formula, δ_c represents the maximum normal gap between the particles when there is non-zero cohesion (m), while α_c represents the maximum tangential gap between the particles when there is non-zero cohesion (m).

When δ_n > δ_c, the cohesion between the particles becomes 0. When the particles are not in actual contact and the distance is less than δ_c, the cohesion reaches its maximum value. The method for calculating the maximum value is shown in Formula (11).

$$F_{pullout}=-\frac{3}{2}\pi ·\gamma ·R^{\ast }$$

(11)

4.3. EDEM Software Simulation Parameter Settings

This study utilized EDEM software for the simulation. The computer had a 13th Gen Intel(R) Core (TM) i5-13490F CPU, 32 GB RAM, and an NVIDIA GeForce RTX 3050 graphics card with 8 GB video memory. The soil model parameters play a crucial role in the simulation results; thus, accurately setting these parameters is of utmost importance. Material parameters, such as the density, Poisson’s ratio, and shear modulus, were determined through the experiments or literature and are listed in

Table 3. Intrinsic parameters of soil and 65 Mn steel.

Material	Density/(kg·m−3)	Poisson’s Ratio	Shear Modulus/Pa
Soil	2.61 × 103	0.41	1.15 × 107
65 Mn	7.865 × 103	0.3	7.90 × 1010

. The contact parameters, including the static friction factor, rolling friction factor, and collision recovery coefficient, were determined based on the data from the measurement tests and relevant references [28,30,31,32]. The recommended ranges for these contact parameters are provided in

Table 4. Factors and levels of the simulation test.

Symbol	Parameter	Low Level	High Level
A	Soil–soil recovery coefficient	0.2	0.8
B	Soil–soil static friction coefficient	0.3	0.9
C	Soil–soil rolling friction coefficient	0.02	0.08
D	Soil–65 Mn recovery coefficient	0.2	0.8
E	Soil–65 Mn static friction coefficient	0.3	0.9
F	Soil–65 Mn rolling friction coefficient	0.1	0.3
G	Surface energy of soil for JKR Model/(J·m−2)	0.2	0.8

. The particle generation method was set to dynamic, with a generation rate of 2000 particles/s and a total of 24,000 generations. The simulation time was set to 18 s, with a time step of 0.2 s and a grid size of 10R. The gravity acceleration was set to 9.81 m/s², and the time step was set to 19% of the Rayleigh time step. The particle size distribution was determined according to Table 3. After 18.0 s, a stable particle pile was formed (Figure 11).

4.4. Plackett–Burman Screening Test and Significance Analysis

The Plackett–Burman test was designed using Design Expert software, with the soil repose angle as the response value to identify the physical parameters that have a significant influence. In the simulation test, eight real parameters (A–G) were designed. Each parameter had two levels, low and high, represented by codes −1 and +1, respectively (Table 4). A center point was also included in the simulation test, resulting in a total of 13 sets of tests [33].

The Plackett–Burman screening test design and results are presented in

Table 5. Design and results of the Plackett–Burman test.

No.	A	B	C	D	E	F	G	Y/°
1	0.8	0.3	0.02	0.2	0.9	0.1	0.8	44.2
2	0.2	0.3	0.02	0.2	0.3	0.1	0.2	30.4
3	0.2	0.3	0.08	0.2	0.9	0.3	0.2	46.7
4	0.2	0.3	0.02	0.8	0.3	0.3	0.8	47.1
5	0.8	0.9	0.02	0.8	0.9	0.3	0.2	33.9
6	0.2	0.9	0.08	0.8	0.3	0.1	0.2	43.9
7	0.8	0.9	0.08	0.2	0.3	0.1	0.8	46.9
8	0.2	0.9	0.02	0.8	0.9	0.1	0.8	47.9
9	0.8	0.3	0.08	0.8	0.9	0.1	0.2	43.3
10	0.2	0.9	0.08	0.2	0.9	0.3	0.8	55.2
11	0.8	0.3	0.08	0.8	0.3	0.3	0.8	48.3
12	0.8	0.9	0.02	0.2	0.3	0.3	0.2	28.9

Note: Y refers to the repose angle.

, with the analysis based on the test results displayed in

Table 6. Analysis of the Plackett–Burman test.

Source	Sum of Squares	df	Mean Square	F-Value	p-Value
Model	674.2	7	96.31	25.46	0.0037 **
A	55.04	1	55.04	14.55	0.0189 *
B	0.9075	1	0.9075	0.2399	0.6499
C	224.47	1	224.47	59.34	0.0015 **
D	12.2	1	12.2	3.23	0.1469
E	55.04	1	55.04	14.55	0.0189 *
F	1.02	1	1.02	0.2699	0.6308
G	325.52	1	325.52	86.06	0.0008 **
Residual	15.13	4	3.78
Cor total	689.33	11

Note: * and ** indicate significance at the 0.05 and 0.01 levels, respectively.

and Figure 12. The order of significance, from largest to smallest, was as follows: surface energy of soil for the JKR model, soil–soil rolling friction coefficient, soil–65 Mn static friction coefficient, and soil–soil recovery coefficient. Consequently, these four influencing factors were chosen for further steepest climbing tests.

Here, p-values less than 0.0500 indicated that the model terms were significant. In this case, A, C, E, and G were significant model terms. Values greater than 0.1000 indicated that the model terms were not significant.

4.5. Steepest Climbing Test Results and Analysis

The range values of A, C, E, and G in Table 5 were too large. To determine the values of each simulation parameter more accurately, a steepest climb test was required. Based on the results of the Plackett–Burman test, the steepest climb test was conducted on the four selected factors that had greater contributions. The relative error between the actual repose angle and the simulated repose angle was calculated to determine the optimal range of the simulation parameters. During the simulation process, the parameters that had little contribution to the repose angle were set to intermediate-level values. Specifically, B was set to 0.6, D was set to 0.5, and F was set to 0.2.

The design and results of the steepest climbing test are presented in

Table 7. Steepest climb test design and results.

No.	G	C	E	A	Repose Angle/(°)	Relative Error/(%)
1	0.20	0.020	0.30	0.80	28.14	21.90
2	0.35	0.035	0.45	0.65	36.97	2.61
3	0.50	0.050	0.60	0.50	47.72	32.45
4	0.65	0.065	0.75	0.35	55.26	53.37
5	0.80	0.080	0.90	0.20	67.31	86.82

. The findings indicate that the repose angle had the smallest relative error in the second set of tests. This suggests that the optimal interval was close to the parameter values used in the second set. Therefore, the parameter values from the first and third groups were considered as the low level and high level, respectively, for the subsequent Box–Behnken tests.

4.6. Box–Behnken Test Results and Analysis

To determine the best parameter combination of G, C, E, and A in the simulation test, a four-factor and three-level test design was conducted based on the results of the climbing test. The repose angle was used as the test index, following the Box–Behnken test principle. A total of 29 tests were performed, and the test plan and results are presented in

Table 8. Box–Behnken test plan and results.

No.	Influencing Factors				Y/°
	G	C	E	A
1	0.5	0.05	0.45	0.65	46.21
2	0.35	0.02	0.45	0.8	31.13
3	0.2	0.04	0.45	0.8	33.51
4	0.35	0.04	0.45	0.65	34.92
5	0.2	0.02	0.45	0.65	31.53
6	0.35	0.04	0.45	0.65	34.88
7	0.35	0.05	0.3	0.65	37.97
8	0.5	0.04	0.45	0.5	36.92
9	0.5	0.04	0.45	0.8	44.94
10	0.35	0.02	0.3	0.65	28.55
11	0.35	0.04	0.6	0.8	40.27
12	0.35	0.04	0.6	0.5	40.94
13	0.35	0.05	0.45	0.5	38.01
14	0.35	0.04	0.3	0.8	33.91
15	0.35	0.05	0.6	0.65	44.86
16	0.35	0.02	0.45	0.5	34.02
17	0.35	0.05	0.45	0.8	42.73
18	0.2	0.04	0.45	0.5	34.18
19	0.5	0.02	0.45	0.65	34.53
20	0.35	0.04	0.45	0.65	34.96
21	0.35	0.04	0.3	0.5	32.32
22	0.5	0.04	0.6	0.65	46.07
23	0.2	0.05	0.45	0.65	37.22
24	0.35	0.04	0.45	0.65	34.94
25	0.2	0.04	0.3	0.65	29.85
26	0.2	0.04	0.6	0.65	34.28
27	0.5	0.04	0.3	0.65	40.73
28	0.35	0.02	0.6	0.65	33.72
29	0.35	0.04	0.45	0.65	34.93

.

Design Expert software was used to perform a multivariate fitting analysis on the test results. A regression model was established between the repose angle and the four independent variables. The regression equation is presented as Equation (12). The results of the regression model variance analysis can be found in

Table 9. ANOVA for the reduced quartic model.

Source of Variance	Sum of Squares	Degrees of Freedom	Mean Square	F-Value	p-Value
Model	657	24	27.4	31,088.73	<0.0001 **
A-G	50.2	1	50.2	57,042.3	<0.0001 **
B-C	60.8	1	60.8	69,047.76	<0.0001 **
C-E	56.1	1	56.1	63,750.11	<0.0001 **
D-A	0.21	1	0.21	240.45	0.0001 **
GC	8.97	1	8.97	10,193.21	<0.0001 **
GE	0.21	1	0.21	235.26	0.0001 **
GA	18.9	1	18.9	21,453.44	<0.0001 **
CE	0.74	1	0.74	840.45	<0.0001 **
CA	14.5	1	14.5	16,452.3	<0.0001 **
EA	1.28	1	1.28	1451.02	<0.0001 **
G2	4	1	4	4541.07	<0.0001 **
C2	0.97	1	0.97	1105.84	<0.0001 **
E2	3.17	1	3.17	3606.98	<0.0001 **
A2	4.78	1	4.78	5434.56	<0.0001 **
G2C	0.4	1	0.4	450.06	<0.0001 **
G2E	3.39	1	3.39	3855.7	<0.0001 **
G2A	5.17	1	5.17	5872.86	<0.0001 **
GC2	0.59	1	0.59	675.06	<0.0001 **
GE2	9.03	1	9.03	10,262.78	<0.0001 **
C2E	1.07	1	1.07	1211.14	<0.0001 **
C2A	0.1	1	0.1	117.63	0.0004 **
CE2	3.09	1	3.09	3508.65	<0.0001 **
G2C2	0.32	1	0.32	369.2	<0.0001 **
G2E2	0.29	1	0.29	334.44	<0.0001 **
Pure error	0.0035	4	0.0009
Cor total	656.60	28
R2 = 1		Adjusted R2 = 1		Adequate precision = 641.18

Note: ** indicate significance at the 0.01 levels.

.

Y = 53.93 − 123.69G − 1587.01C + 118.71E − 45.53A + 6390.00GC + 388.70GE − 236.85GA − 1608.15CE + 373.70CA − 25.11EA + 234.89G² + 26932.96C² − 179.36E² + 47.37A² − 6562.96G²C − 1350.37G²E + 476.30G²A − 94962.96GC² − 120.49GE² − 21629.63C²E + 6740.74C²A + 3681.48CE² + 112592.59G²C² + 1071.60G²E²

(12)

The regression model shows a significant result with a p-value of less than 0.0001. The coefficient of determination R² was 1, the adjusted R² was 1, and the adequate precision was 641.18, indicating the high accuracy of the model. This suggests that the table effectively illustrates the relationship between the repose angle and the four factors.

4.7. Regression Model Interaction Effect Analysis

In this experiment, the soil repose angle was used as the response variable. The data were analyzed using Design Expert software through multiple regression fitting. The response surface, as depicted in Figure 13, was generated to visualize the relationship between the influencing factors and the response value.

Figure 13 illustrates that the influence of G, C, E, and A on the repose angle aligns with the findings of the Plackett–Burman test mentioned earlier. Notably, the response surface curve for G exhibited the sharpest incline, suggesting it had the most pronounced effect on the repose angle. Conversely, the curve for A showed a more gradual slope, indicating that its effect on the repose angle was the least significant. Furthermore, G, C, and E had a positive effect on the repose angle, whereas A had a negative effect.

4.8. Parameter Optimization and Simulation Verification

In the optimization module of Design Expert 11.0 software, the regression equation needs to be optimized and solved with a target value of 36.03°. The optimization target value and constraints are represented in Equation (13).

$$\{\begin{array}{c}36.03\\ s.t.\{\begin{array}{c}0.2\le \mathrm{G}\le 0.5\\ 0.02\le \mathrm{C}\le 0.05\\ 0.3\le \mathrm{E}\le 0.6\\ 0.5\le \mathrm{A}\le 0.8\end{array}\end{array}$$

(13)

The optimal parameter combination obtained through the optimization solution consisted of a value of G at 0.400, a value of C at 0.040, a value of E at 0.404, and a value of A at 0.522. To verify the reliability and authenticity of the discrete element simulation parameters post-calibration of the soil, a simulation test was conducted using the optimal parameter combination as the EDEM simulation parameters. The repose angle of the soil was measured at 36.91° with a relative error of 2.44%. For verification purposes, a comparison of the test results is presented in Figure 14 to ensure the reliability and authenticity of the simulation test.

The reliability of the calibration parameters was further confirmed using the 65 Mn steel rod insertion method based on the optimal parameter combination identified previously (Figure 15). The dimensions of the container in which the soil was placed (base area and height) corresponded to the dimensions of the glass container. In both the physical and simulation tests, a 3 mm wide and 100 mm long steel rod was dropped freely from a height of 150 mm onto the soil, and the depth of the steel needle inserted into the soil after the fall was measured. The physical test and simulation test were repeated eight times each, and the statistical results are presented in Figure 16. The final result was determined by taking the average value from both the physical and simulation tests. The physical test yielded an average insertion depth of 19.31 mm, while the simulation test yielded an average insertion depth of 18.98 mm, resulting in a relative error of 1.71%. This test once again verified the reliability of the calibrated parameters.

5. Discussion

In numerous studies, the discrete element calibration method for discrete state materials primarily relies on the repose angle as a response variable. While the calibration process remains largely consistent across studies, the methods of verification vary. Some researchers have opted to use the physical repose angle under optimal conditions for verification, while others have utilized the resistance of material-contacting components or the movement of materials in the device [34,35,36]. Our study employed verification through the physical angle of repose measurements and steel rod drop tests. It is important to note that calibration results may differ based on variations in soil quality, environmental conditions, particle size, and initial parameter ranges. According to the relevant theories, the material–material static friction coefficient, material–material rolling friction coefficient, material–contact component static friction coefficient, material–contact component rolling friction coefficient, and surface energy of soil for the JKR model are all positively correlated with the repose angle. On the other hand, the material–material collision recovery coefficient and material–contact component collision recovery coefficient are inversely related to the angle of repose.

The measurement of the repose angle lacks a standardized method, resulting in inaccuracies that can affect the significance of the influencing factors. Even factors expected to positively impact the repose angle may be distorted by bias, yielding misleading results. For example, in the study by Deli J et al. [37], the static friction coefficient of the cotton stalk–cotton stalk had a negative impact on the repose angle (it should have been positive). In our study, due to the existence of deviation, the soil–soil static friction coefficient in Figure 12 showed a negative effect.

6. Conclusions

The Hertz–Mindlin with JKR contact model was utilized to calibrate the discrete element parameters of Camellia oleifera forest soil in Changsha, Hunan, China, based on discrete element EDEM simulation software. By combining the physical tests and simulation tests, along with Design Expert software, the Plackett–Burman test identified the key factors affecting the soil repose angle, which were the surface energy of the soil for the JKR model, soil–soil rolling friction coefficient, soil–65 Mn static friction coefficient, and soil–soil recovery coefficient. A regression model was established through the Box–Behnken test to analyze the variance and interaction effects of the four factors on the repose angle. The optimal parameter combination, determined to be the surface energy of the soil for the JKR model (0.400), soil–soil rolling friction coefficient (0.040), soil–65 Mn static friction coefficient (0.404), and soil–soil recovery coefficient (0.522), was validated through stacking simulation tests. The relative error between the optimal parameter combination and the actual physical repose angle was found to be 2.44%, confirming the reliability of the soil discrete element parameters. Furthermore, validation through a steel rod insertion test showed a relative error of 1.71% between the insertion depth in the simulation test and the physical test, further supporting the reliability of the soil discrete element parameters.