DEM-based calibration of paddy soil contact parameters for improved soil-tool interactions

1 Introduction

Rice is one of Sierra Leone's staple crops, and paddy cultivation is a critical component of rice production (Tienan et al. 2023). Precise paddy field preparation procedures improve rice yield, irrigation uniformity, water and fertilizer use efficiency, and control weeds, insects, and diseases, promoting sustainable agricultural growth (Asmarani et al. 2023; Xing et al. 2024; Korba et al. 2024). Understanding the interaction between agricultural implements and paddy soils is key to optimizing agricultural implements and minimizing soil disturbance in paddy fields. Analytical, experimental, and numerical methods have been used to study the interaction between soil and tools (Ucgul 2023; Khosravani et al. 2023; Hashaam et al. 2023). Analytical methods are limited to simple geometries and assumptions about soil failure, while experimental methods (including field tests and laboratory soil bin tests) are time-consuming and expensive (Wang et al. 2022a; Aikins et al. 2023). Numerical methods can be used to investigate complex tool geometries and a wider range of soil dynamic parameters (Wang et al. 2022b). Discrete element modelling (DEM) can be used to simulate the mechanical behavior of discontinuous particles, such as soil, without restrictions on the size of particle displacement or tool design (Du et al. 2022; Zhang et al. 2022; Qiu et al. 2023). The discrete element method (DEM) has become a popular modelling approach in agricultural research in recent years. In 2022, Wang et al. used DEM to create a deep tillage tool-soil interaction model that accurately predicted the microscopic behavior of soil at various locations. Dilip et al. (2022) used DEM to study the effect of tine spacing on soil disturbance behavior; they found that a tine spacing of 400 mm performed better than other spacings. Kim et al. (2022) used DEM simulations to investigate how the operating depth of deep tillage tools affects their performance in clay soils. Zhou et al. (2022) used the discrete element model (DEM) to establish a paddy soil model with slump expansion as the performance indicator. They found that soil surface energy, soil-soil rolling friction coefficient, and soil-steel static friction coefficient significantly affected the overall slump expansion error between simulation and verification results. Wang et al. (2023) used DEM to simulate soil sweeping interactions and conducted field verification. Song et al. (2025) used DEM to simulate the soil loosening process induced by vibratory deep loosening equipment. Previous DEM research has primarily focused on simulating soil disturbance behavior at relatively low water content (25%). In contrast, paddy soils have significantly higher water content (>40%), viscosity, and fluidity, resulting in significantly different disturbance behaviours. However, previous literature has provided fewer DEM parameters, which are crucial for simulating paddy soil-tool interactions (Chen et al. 2023).

This study expands upon existing DEM literature, aiming to provide a more accurate simulation model for studying the interaction mechanisms between agricultural implements and paddy soil. Paddy soil slump was used as the experimental metric, based on slump tests. Test factors included Poisson's ratio, shear modulus, soil surface energy, soil-soil restitution coefficient, soil-soil static friction coefficient, soil-soil rolling friction coefficient, soil-steel restitution coefficient, soil-steel static friction coefficient, and soil-steel rolling friction coefficient. These DEM parameters were calibrated using the Plackett-Burman test, steepest ascent test, and Box-Behnken test in Design-Expert 13. Comparing the slump test and simulated values helped assess the accuracy of the DEM parameters during the calibration process. The results of this study provide a basis for selecting ideal DEM parameters for paddy soil models, thereby facilitating future research on the relationship between soil and agricultural implements to provide optimal working depth and minimize soil disturbance.

2 Materials and methods

2.1 Soil slump experiment

This study used the slump test as a basic test for calibrating paddy soil parameters. The paddy soil type is sandy loam with an average water content of 35.66%. The slump was used to determine the physical parameters of the paddy soil, and the slump error was used as a test indicator to calibrate the paddy soil model parameters. A slump cone was used to measure the slump height of the paddy soil. The cone was supported by a horizontal base plate and had a height of 300 mm, with upper and lower diameters of 100 mm and 200 mm, respectively. After loading, the paddy soil was compacted, and excess soil was scraped off the top of the cone. The cone was then raised at a constant speed over a period of 5 seconds. After the soil slumped, the height of the cone (h) was measured (Figure 1), and the average of five repeated measurements was recorded.

2.2 DEM contact model

Discrete element modeling (DEM) software was used to run the simulation. The Hertz-Mindlin and Johnson-Kendall-Roberts (JKR) models are commonly used for analyzing wet cohesive soils (Zhou et al. 2022; Tienan et al. 2023). This model is based on the Hertz-Mindlin (no-slip) contact model, which considers the effects of inter-particle plastic deformation and inter-particle adhesion on particle motion (Hensh et al. 2022; Jing et al. 2023; Xu et al. 2023; Siyuan and Duruo 2023). Jing et al. (2023) effectively used this contact model to simulate cohesive soils. Because paddy loam is extremely cohesive, the Hertz-Mindlin and JKR models were selected as the contact models for the discrete element soil model in this study. In the JKR contact model, the normal elastic contact force FJKR between particles is determined as follows:

The relationship between the contact radius between particles r and the degree of overlap δ, is determined as:

The equivalent contact radius R' and elastic modulus E' are determined as follows:

where, F_JKR is the JKR normal contact force, N; γ - soil surface energy, J.m^-2; E ' - equivalent modulus of elasticity; R ' - the equivalent contact radius, mm; V₁ and V₂ are the particle Poisson’s ratios; R₁ and R₂ are the particle radii in mm.

2.3 Model parameter calibration

Material and interaction properties are the two basic categories of soil parameters in DEM. Material properties are the basic characteristics of the soil, such as bulk density, shear modulus, particle size and shape, and Poisson's ratio. The characteristics a particle exhibits when it comes into contact with boundaries, surfaces, and other particles are known as interaction properties (Wang et al. 2023; Yuan et al. 2023). Such properties include surface energy, coefficients of restitution, static and rolling friction. According to a combination of initial testing and the literature, the nominal soil particle radius was 1 mm (Zhou et al. 2022; Aikins et al. 2023). The bulk density of the soil was obtained as 1400 kg.m^-3. According to earlier studies the shear modulus, Poisson's ratio, and density of steel were calculated. The soil's shear modulus, surface energy, and Poisson's ratio, as well as the soil-soil and soil-steel coefficients of restitution, static, and rolling friction, were chosen from the literature (Table 1) (Song et al. 2022; Zhou et al. 2022; Wang et al. 2022; Tienan et al. 2023; Aikins et al. 2023).

2.4 Slump simulation test

Based on the actual dimensions of the slump cone used in the slump test (Figure 1), a 1:1 scale 3D model was constructed and imported into the DEM. The slump test procedure was followed when simulating paddy soil, and the slump height of the paddy soil was measured using the post-processing features of the DEM. Using the slump height error as the test metric, the Plackett-Burman test was used to identify simulation factors that significantly affected the slump error. The steepest ascent test was used to further narrow the range of optimal values for these key factors. The Box-Behnken test was then used to develop a regression model between the test parameters and the test metric (Yuan et al., 2023). To find the ideal parameter combination, the regression equation was optimized and the results were verified by simulation. A screenshot of the particle stacking state during the slump simulation is shown in Figure 2.

The slump height error ( δ r ) was employed to determine the variance between the simulated and test slumping values. The slump error is calculated as follows:

where, h_e - represent test slump value, mm; h_s - represent simulated slump value, mm

3 Results and discussions

3.1 Slump test results

The slump test was conducted five times, with slump values ranging from 210.77 mm to 229.5 mm, and an average of 223.92 mm, as shown in Table 2. The coefficient of variation of the slump values was 3.39%, indicating that the test results were highly accurate.

“C_V” in the table represents coefficient of variation

3.2 Plackett-Burman test

When there are many experimental factors, the Plackett-Burman test can quickly identify those with significant impact on an indicator. Nine simulation parameters were selected for this experiment, with each parameter taking values at high and low levels (+1 and -1), as shown in Table 1. Using Design-Expert software, a Plackett-Burman test design was performed according to the parameter levels in Table 1. Table 3 shows the experimental plan and results. An analysis of variance was performed on the experimental data. Table 4 highlights the significance of each parameter's impact on the test indicator (slump height, h). The nine experimental factors are arranged in descending order of contribution to the test indicator. According to Table 4, x₂, x₆, and x₈ contribute relatively significantly to the total relative error (δ_r). Therefore, the steepest ascent test and Box-Behnken test were used to calibrate these three parameters. The p-values for x₂, x₆, and x₈ were all less than 0.05, indicating that they have an impact on the test indicators. Zhou et al. (2022) also found that, despite using a high paddy soil moisture content of >40%, soil surface energy, soil-soil rolling friction coefficient, and soil-steel static friction coefficient had a significant impact on soil slump error.

3.3 Steepest ascent test

To further investigate the ranges of soil parameters of x₂, x₆, and x₈, steepest ascent test was designed and conducted to examine their impact on the total relative error (Table 4). This test was used to screen and investigate the effective ranges of these three key factors. Table 5 shows that the slump height error decreases with increasing levels of each factor. The minimum relative error (6.98%) occurs at the fourth simulation level, indicating that the ranges of x₂, x₆, and x₈, are relatively reasonable: 0.6–1.10, 0.055–0.11, and 0.5–0.9, respectively. With the fourth simulation level as the center point, the third and fifth simulation levels serve as parameter intervals for further calibration. The resulting feasible range of soil parameters was used in subsequent Box-Behnken test.

3.4 Box-Behnken test

The Box-Behnken test method was used to examine the specific effects of soil surface energy (x₂), soil-soil rolling friction coefficient (x₆), and soil-steel static friction coefficient (x₈) on the slump value (h), as well as the relative error between the simulated and measured h values. The test determined the factor set that optimizes the test indicator, the slump error of paddy soil. Table 6 shows the coding of the test factors.

The Box-Behnken test was conducted using the parameters listed in Table 6. Table 7 shows the test plan and results. Multiple regression analysis was performed on the test results using Design-Expert software to obtain a slump error regression model and perform analysis of variance (ANOVA).

Table 8 presents the ANOVA results of Box-Behnken test. The results showed that the effect of x₆ on the slump error is not significant (P > 0.05), while the effects of x₂, x₈, x₂x₆, x₂x₈, x₆x₈, x₂², x₆², and x₈² are all very significant (P < 0.01). The P values of the regression model and the lack of fit term of the slump error fitting are less than 0.01 and greater than 0.05, respectively, indicating that the model fits well. The coefficient of determination R² of the regression equation was 0.9885, and the adjusted coefficient of determination adj-R² was 0.9706. The results indicated that the regression model can be used to predict the slump of paddy soil. The optimized slump error regression equation was calculated as follows:

3.5 Parameter optimization and verification

The regression model was optimized using the Design-Expert optimization tool, taking the range of test factor values as boundary conditions and the slump error as the target parameter to determine the optimal parameter combination affecting the test index. The optimization results showed a soil surface energy of 0.869 J·m⁻², a soil-soil rolling friction coefficient of 0.109, and a soil-steel static friction coefficient of 0.651. Five slump tests were conducted using the optimized x₂, x₆, and x₈ combinations in discrete element simulations to verify the optimization results. Simulations were conducted to validate the optimized results for the ideal parameter combination (Table 9). The average error between the simulated and experimental values for slump height was 5.59%, with a coefficient of variation of 3.49%. Zhou et al. (2022) and Tienan et al. (2023), despite using different paddy soil properties, also found similar coefficients of variation. These results suggest that the paddy soil discrete element model using this parameter combination can provide fundamental information for simulation studies of soil-tool interactions.

Conclusion

This study used a combination of the Plackett-Burman test, the steepest ascent test, and the Box-Behnken test to determine the DEM parameters for paddy soil using the JKR DEM model. The accuracy of the DEM model was evaluated using a slump test. For a given paddy soil model, soil surface energy (x₂), soil-soil rolling friction coefficient (x₆), and soil-steel static friction coefficient (x₈) significantly influenced the model accuracy. Given an overall relative error of 5.59% between the slump test and the simulated values, the ideal combination of x₂ of 0.869 J.m^-2, x₆ of 0.109, and x₈ of 0.651 was obtained. These results demonstrated the accuracy of the DEM model, which can be used to establish paddy field soil simulation models and design high-performance soil tools.