Modelling of soil-tool interaction allows to gain an in-depth understanding of tool performance and optimise soil engaging tools without performing expensive and time-consuming field tests (Fielke, 1999; Shmulevich et al., 2007; Ucgul et al., 2014a). The discrete element method (DEM) is a dis-continuum numerical method to model the mechanical behaviour of discontinuous soil particles. DEM simulations can be performed in two and three-dimensions (Asaf et al., 2007; Shmulevich, 2010). For 2D DEM, real soil particle sizes can be used and thus the accuracy of the simulation results can be high (Ucgul et al., 2014b), for instance the modelling of the soil loosening process caused by a vibrating subsoiler by Tanaka et al. (2000); the modelling of the cutting blades by Zhang & Li (2006) and Zhang et al. (2008); and the modelling of a pendulum type cutting blade test by Momozu et al. (2003). However, in 3D DEM, using real soil particle sizes is not practical due to much larger number of the particles required to assembly a soil bin and the computational limitations. To gain solutions in a timely manner, many 3D DEM studies selected particles larger than the particle sizes found in soil (Chen et al., 2013; Ucgul et al., 2014a; Hang et al., 2018; Li et al., 2018; Tamas, 2018).

To date, a few attempts have been made to examine the effect of soil particle size on soil disturbance behaviors and cutting forces, e.g. appropriate particle radii gave a higher level of agreement between simulated and measured soil disturbance profiles (from a sweep tool impact) (Milkevych et al., 2018); reducing particle radii from 10 mm to 5 mm gave more accurate predictions of cohesion and friction angle between particles (Ucgul et al., 2015); simulations with a 20 mm thick layer of 1.5 mm radii particles over a base of 10 mm radii particles provided a more accurate furrow profile (Ucgul et al., 2014a); the mean particle size had a negligible effect on the small strain shear modulus for mono-sized soils (Gu et al., 2017); the magnitude of internal friction angle increased as the particle size distribution was made wider (Coetzee & Els, 2009). Previous DEM works in relation to soil tillage (Ucgul et al., 2014a, 2015) mainly focus on the effect of soil particle size under a particle-size distribution (PSD) which was mostly selected arbitrarily. Although good correlations have been obtained between average simulated and measured tillage forces, some important soil disturbance characteristics or soil cutting forces in real time have either not been provided or not well correlated with the measured results.

The significance and mechanism for the effect of PSD on soil disturbance characteristics and cutting forces were not well documented in previous DEM tillage studies. Therefore, the objectives of this study were to (1) develop DEM models with different PSDs using EDEM 2.7 software (DEM-Solutions Inc., UK), (2) investigate the effect of PSD on predicted soil disturbance characteristics and cutting forces, and (3) validate the simulation results using the laboratory soil bin data.

Soil bin studies

Description of the equipment and the soil bin

The subsoiler tested (Fig. 1a) was 159-mm wide. It has a pair of wings and an arc-shaped shank selected based on the Chinese standards (JB/T 9788-1999) (Hang et al., 2018). It was mounted on a toolbar at a rake angle of 23° (Fig. 1). The soil bin tests were conducted at Northwest A&F University, Yangling, China. The study area is characterized by a temperate semi-arid continental monsoon climate. The soil bin contained a Lou soil (15.67±0.75% gravel, 74.61±1.75% sand, 9.13±1.63% silt and 0.58±0.09% clay) that developed on parent loess (Huang et al., 2016; Hang et al., 2017). According to the physical parameters of the selected field soil (e.g. density and moisture content), the soil for the soil bin test was prepared utilizing a layered method. Initially, the top layer soil (170 mm from the soil surface) was removed carefully. The rest soil was watered, then a vibratory rammer (frequency of ramming: 7-11 Hz) and a steel roller were used to compact the soil in sequence. After preparing the hardpan soil, the removed soil was evenly backfilled. And then similar operations, i.e. spraying water and compaction, were conducted to prepare the top layer soil.

FigureFigure 1.The subsoiler (a) and soil bin cart (b) with the three-point hitch.

Measurements

After soil preparation, soil cores were taken from undisturbed areas in the bins for three repetitions. Soil samples were brought to the laboratory to measure the soil moisture content and bulk density. The thickness, soil moisture content and bulk density of the top layer were 0.17 m, 17% and 1400 kg m-3, respectively. The corresponding values for the hardpan were 0.13 m, 22% and 1830 kg m-3. The subsoiler was run in the soil bin at a constant speed of 3 km h-1 and a working depth of 300 mm in three replicates (Fig. 1b). Draft and vertical forces of the subsoiler were measured by sensors installed between the soil bin carriage and the three-point hitch of the toolbar.

Shear and compression forces exerted on the soil by a tool is the major factor that causes the soil structural failure (Perfect et al., 2002; Hang et al., 2018). Soil rupture distance ratio is a critical parameter to characterise the soil structural failure (Hettiaratchi et al., 1966; Godwin & Spoor, 1977). A tillage tool with smaller soil rupture distance ratio has lower tractor power requirement due to smaller soil cutting forces. Soil rupture distance ratio was therefore often used to evaluate the working performance of soil engaging tools, e.g. the study of Yang et al. (2018). As per Hettiaratchi et al. (1966) and Yang et al. (2018), soil rupture distance ratio was calculated as,

where, m is soil rupture distance ratio; f is soil rupture distance (Fig. 2b); d is working depth of a tool.

Smaller height of accumulated soil is often desired for an ideal subsoiler as it corresponds to a better sliding-cutting performance of a subsoiler (Zhou et al., 2019). The height of accumulated soil (hs) was defined by the height of disturbed soil on the original surface before the shank of the subsoiler (Fig. 2b). The soil rupture distance ratio and height of accumulated soil were measured using an I-SPEED TR high-speed video camera (Olympus Co., Japan) (Fig. 2a). The vertical distances between various points on the shank and the lowest position of the subsoiler were f irstly marked before the tests (Fig. 2c). During tillage, the high-speed video camera was used to record the time interval between the moment when the soil on the white line (vertical to the center of subsoiler path) began to move (i.e. Moment 1) and the moment when the subsoiler arrived at the first of the three white lines (i.e. Moment 2) (see Fig. 2b); and then soil rupture distance (f) was calculated from the travel speed (3 km h-1) and the time interval. The working depth is the distance between soil surface and the lowest position of the furrow which was excavated manually after tillage. The height of accumulated soil was determined by the vertical distance captured during tillage (Fig. 2d) and the working depth.

There is an increasing concern about soil compaction and hardpan formation with the successive use of shallow tillage tools and larger tractors (Chen et al., 2005; Zeng et al., 2017). Hardpan disruption helps to restore soil bulk density and improve water infiltration, root penetration, and crop growth (Shahgoli et al., 2009; Zeng et al., 2017). Soil density change rate of the hardpan was used to evaluate the subsoiler’s performance in terms of the quality of hardpan loosening. It was calculated as follows:

where, P is soil density change rate, %; ρb and ρa are densities of hardpan soil before and after tillage, respectively, kg m-3.

Measurements were performed at three random locations for soil bulk density and six random locations for working depth in each plot. The experiment was repeated three times. In total, three groups of draft and vertical forces, soil rupture distance, height of accumulated soil and soil bulk density, and eighteen working depths were collected and averaged for the later validation of DEM simulations.

FigureFigure 2.The high-speed video camera (a) for measuring soil rupture distance and height of accumulated soil (the black arrow stands for the direction of travel and the centre of subsoiler path); diagram (b) showing soil rupture distance f and height of accumulated soil hs during tillage (d stands for working depth of the subsoiler; A and E stand for intersections of the first of the three white lines and the centre of subsoiler path at moment 1 and moment 2, respectively); subsoiler (c) with a marked shank; and the moment (d) when the subsoiler arrived at the first of the three white lines.

DEM simulations

DEM contact model and parameters

Hertz-Mindlin with bonding (HMB) model was employed in this study. In the HMB model, a “cylindrical cementitious material” (i.e. bond) was added to represent the liquid bridge between particles (Fig. 3). The bond can withstand forces and moments whose magnitudes depend on the micro-properties (e.g. critical normal and shear stresses) of the bond (Chen et al., 2013; Ding et al., 2017). In this study, the bond stiffness was 5e7 N m-3; the bond critical stresses of top layer and the hardpan were 30 and 40 kPa respectively based on the publish data (Wang et al., 2018). The bond radius was determined in accordance with the particle nominal radius, soil density and moisture content (Ding et al., 2017; Wang et al., 2018). The DEM parameters can be grouped into two categories, namely material and interaction properties. The material properties were obtained from a combination of measurements and published data (using similar soil conditions). The soil density was obtained by measurement (using a precision scale). The shear modulus and Poisson's ratio of soil and the density and shear modulus of steel used in this study were published data (Huang et al., 2016; Wang et al., 2018), as shown in Table 1. The interaction properties mainly consisted of the coefficient of restitution, the coefficient of rolling friction and the coefficient of friction between materials. The coefficient of rolling friction between soil particles and between soil particles and the steel were obtained by calibration based on the angle of repose test and the inclined plane test, respectively. The coefficient of restitution (0.6) and the coefficient of friction between materials were obtained from previous studies (Ding et al., 2017; Shi et al., 2017; Hang et al., 2018; Wang et al., 2018). A typical percentage of Rayleigh time step of 20% (Liu et al., 2016) was selected in the simulations. The time step was automatically calculated by EDEM software according to particle sizes and the given percentage of the Rayleigh time step.

FigureFigure 3.The soil-subsoiler interaction model.

Soil-subsoiler interaction model

To avoid the effect of bin walls on soil particle flows during subsoiling, the vertical and horizontal dimensions of the virtual bin were set larger than tillage depth and soil disturbance width, respectively. The longitudinal dimension of the bin was determined in accordance with the computation power of the computer and subsoiler length; i.e., the bin should allow for the subsoiler to reach a stable condition in terms of the draft force (Mak et al., 2012). According to above criteria, dimensions of the virtual soil bin were set as 1.0× 0.6 × 0.4 m (length × width × depth). For the soil bin, the depth of 0-0.17, 0.17-0.3, and 0.3-0.4 m below soil surface were top layer, hardpan, and subsoil, respectively (Fig. 3). The subsoiler model was developed using CATIA V5R20 software, and then it was positioned at one end of the soil bin at a constant working depth of 300 mm prior to simulations, as shown in Fig. 3.

In HMB model, soil particle shapes are less important since soil particles are bonded together and move as ‘‘ag-gregates’’ (Chen et al., 2013). Considering the aggregate size ranging from 1 to 49 mm for both fine and coarse soils (Mak et al., 2012), many researchers used 8 mm or larger radii spheres as the soil particle models to gain DEM simulations in a timely manner (Zheng et al., 2016; Ding et al., 2017; Ucgul et al., 2017; Wang et al., 2018). To investigate significance and mechanism for the effect of PSD on soil-subsoiler interactions in DEM models, si-mulations were run using spherical particles with nominal radii ranging from 5 to 15 mm with an interval of 2 mm. To describe the particle size range concisely, PSD was de-fined as follows:

where Rn is particle nominal radius; Rmax and Rmin are maximum and minimum soil particle radii in the virtual soil bin, respectively.

For a given Rn, the larger the PSD, the wider the particle size range. PSD in most previous DEM simulations ranged from 0 to 0.4 (Shmulevich, 2010; Tamas et al., 2013; Ucgul et al., 2014a, b; Li et al., 2016; Zeng et al., 2017; Yang et al., 2018). Based on the previous DEM studies and compu-tation power of the computer used for this study, the range of PSD tested was determined: 0, 0.1, 0.2, 0.3, 0.4, 0.8, and 1.2. Each simulation was repeated three times and a total of 126 simulations were performed in 42 virtual soil bins.

Table 1.Major DEM model parameters.

Data collection and analysis

Draft and vertical forces from each simulation were mo-nitored over the entire travel of the subsoiler using EDEM software. Before the subsoiler entered the virtual soil bin, the impact from the subsoiler on particle velocities is a progressive process. When the subsoiler was in the cen-tre of the soil bin, the longitudinal-sectional view of the bin was obtained using the “Clipping” module of EDEM software. The soil rupture distance of subsoiler on the soil surface (f) and the height of accumulated soil (hs) were obtained through the instantaneous particle velocity field shown in Fig. 4; and then soil rupture distance ratio (m) was calculated.

A measuring box was placed in the centre of the soil bin and within the hardpan prior to the simulation (Fig. 5); the mass of particles in the box before and after tillage was then collected. Soil density change rate of the hardpan (P) was calculated using the collected mass of particles according to Eq. (2) and Eq. (5).

where, P is soil density change rate, %; mb and ma are the mass of particles in the measuring box before and after tillage, respectively; V is the volume of the measuring box.

For a given particle nominal radius examined, analysis of variance (ANOVA) was conducted using SPSS 19 statistical software (IBM Corporation, Armonk, NY, USA) to examine the effects of PSD (experimental factor) on soil rupture distance ratio, height of accumulated soil, soil density change rate, and soil cutting forces (draft and vertical). Means of variables were compared between PSDs using Duncan’s multiple range tests to detect differences of the variables between any two treatments. The probability level for the analyses was 5%. The method of statistics can be used to investigate the significance for the effect of PSD on soil disturbance characteristics and cutting forces, and guide the decision of particle-size distribution in the soil-tool interaction DEM models.

FigureFigure 4.Particle velocity field showing the soil rupture distance (f) and the height of accumulated soil (hs).

FigureFigure 5.Particle velocity field showing the soil rupture distance (f) and the height of accumulated soil (hs).

Soil rupture distance ratio

Figure 6 shows the variation of soil rupture distance ratio (m) with the increase of PSD from 0 to 1.2 when different nominal radii (Rn) are used in the discrete element models. Effects of PSD on the m value were significant when particles of 9 mm nominal radii or larger were used in the discrete element model (p < 0.05) and no particular trend was found with the increase of PSD from 0 to 1.2 (Fig. 6). However, with Rn = 7 mm or smaller, there was no statistical differences detected in soil rupture distance ratio, which remained fairly constant over the entire range of PSD examined. For the models with particles of 9 mm nominal radii or smaller, with the increase of PSD, the predicted m had lower relative errors (< 11.6%) as compared with the experimental data measured using the high-speed video camera. For the models with particles of 11 mm nominal radii or larger, with the variation of PSD, the lowest relative error for predicted m was less than 1.9% and the highest relative error was larger than 17.3%. For the models with particles of 7 mm nominal radii or smaller, the simulated soil rupture distance (f) was in close agreement with that measured in the laboratory (422.9 mm) over the entire range of PSD examined (Fig. 7). However, with the variation of PSD, the predicted soil rupture distance varied in a much wider range for the models with particles of greater nominal radii (>7 mm).

Height of accumulated soil

For all nominal particle radii examined, the height of accumulated soil (hs) varied in a wide range when particle-size distribution (PSD) was increased from 0 to 1.2 (Fig. 8). For models with particles of 9 mm nominal radii or smaller, there were no statistical differences found in hs among the tested PSDs (p> 0.05) (see Table 2). However, with further increase in Rn (11 or 13 mm), effects of the PSD on hs were significant (p< 0.05). For models with particles of any given nominal radius examined, with the variation of PSD, the lowest relative error for the predicted hs was less than 6.9% and the highest relative error was larger than 13.1%, as compared with the measured results.

Soil density change rate

Figure 9 demonstrates how the soil density change rate (P) varies during a simulation run of the subsoiler. At time t1, the subsoiler began to move and no variation of the P was observed as the subsoiler did not contact with particles in the measuring box. At time t2, the P value increased quickly as the subsoiler’s cutting share began to enter the measuring box. With further increase in time (t3), the P value increased to a peak when the subsoiler was fully in contact with the measured particles and the hardpan had the maximum soil porosity. Finally, the P value dropped to a constant value as the soil settled down. With Rn=11 mm or smaller, increasing particle-size distribution (PSD) from 0 to 1.2 did not make any differences in soil density change rate (P) (p> 0.05) (Fig. 10). However, for the models with par-ticles of 13 mm nominal radii or larger, the effect of PSD on the P was significant (p < 0.05). With Rn = 5 mm, the predicted P did not vary much with the increase of PSD and had low relative errors (< 7.4%) as compared with the measured value. For the model with particles of any given nominal radius of 7 mm or larger, with the variation of PSD, the lowest relative error for predicted P was less than 8.8% and the highest relative error was larger than 21.3%.

FigureFigure 6.Soil rupture distance ratio (m) as affected by particle-size distribution (PSD) when different particle nominal radii are used in the virtual soil bin. R0: laboratory soil bin data; R5-R15: particle nominal radii (Rn) ranging from 5 to 15 mm; different letters (a-c) mean significant difference at p < 0.05 for a given particle nominal radius; error bars represent standard deviations from the replicates.

FigureFigure 7.Velocity fields of the longitudinal soil failure as affected by particle-size distribution (PSD) when different particle nominal radii (Rn) are used in the virtual soil bin: a) Rn = 5 mm, b) Rn = 7 mm, c) Rn = 9 mm, d) Rn=11mm, e) Rn = 13 mm and f) Rn = 15mm

FigureFigure 8.Height of accumulated soil (hs) as affected by particle-size distribution (PSD) when different particle nominal radii are used in the virtual soil bin. R0: laboratory soil bin data; R5-R15: particle nominal radii ranging from 5 to 15 mm; different letters (a-c) mean significant difference at p < 0.05 for a given particle nominal radius; error bars represent standard deviations from the replicates

Table 2.Significance (p value) for the effect of PSD on different variables for a given particle nomi-nal radius (Rn) based on the ANOVA outputs.

FigureFigure 9.An example of soil density change rate (P) from a model with particle nominal radii of 7 mm and particle-size dis-tribution of 0.1.

Draft and vertical forces

The draft and vertical forces of the subsoiler were calculated from the steady portion of the force curves (Fig. 11). At time t1, both draft and vertical forces were gradually increased. With the increase of time, both draft and vertical forces remained fairly constant when the subsoiler was in the middle of the soil bin (t2). Ulti-mately a large spike was produced as the subsoiler was going through the wall (t3). Figure 12 demonstrates how the particle-size distribution (PSD) affects the predicted draft and vertical forces. Effects of PSD on the predicted vertical and draft forces were significant (p< 0.05) (see Table 2). Vertical and draft forces had different trends on aspects of their variations with the PSD. For the model with particles of any given nominal radius examined, an overall increase trend was detected for draft force with the increase of PSD from 0 to 1.2. In contrast, vertical force fluctuated around a constant value and no trend was detected. For the nominal particle radii: 5, 7, 9, 11, 13, and 15 mm, the predicted average vertical forces of all PSDs examined were -907.95, -816.83, -617.46, -320.93, -146.49, and 22.52 N, respectively; i.e., increasing particle radii generally gave lower vertical forces.

For the model with particles of any given nominal radius examined, with the variation of PSD, the lowest relative error for predicted draft force was less than 8.5% as compared with the measured results; however, the lowest relative error for predicted vertical force was small (<25%) only when small nominal radii (<9mm) particles were used in the models, regardless of the magnitude of PSD. Above results implied that for the given soil, more accurate draft force can be achieved by selecting an appropriate PSD in DEM models; however, the predicted vertical force can correlate well with experimental results only when particles of smaller nominal radii (7 mm or smaller) are used in the discrete element model.

FigureFigure 10.ASoil density change rate (P) as affected by particle-size distribution (PSD) when different particle nominal radii are used in the virtual soil bin. R0: laboratory soil bin data; R5-R15: particle nominal radii ranging from 5 to 15 mm; different letters (a-c) mean significant difference at p < 0.05 for a given particle nominal radius; error bars represent standard deviations from the replicates

FigureFigure 11.Typical simulated draft and vertical forces monitored by EDEM software.

Accurate predictions of soil disturbance behaviors and cutting forces are seriously essential, which enable researchers and engineers to gain insight into the performance of soil engaging tools. For instance, accurate soil rupture distance ratio (m) is critical for the force prediction and characterizing the soil structural failure (Godwin & Spoor, 1977; Yang et al., 2018); height of accumulated soil (hs) is an important parameter to characterise the soil furrow profile which reflects the longitudinal accumulation degree of soil on the surface after the subsoiling operation (Huang et al., 2016); soil density change rate (P) can be used to evaluate the tool performance on aspect of the quality of soil loosening (Zheng et al., 2016). The draft and vertical forces are two important parameters for any soil-engaging tool, as they determine the tractor power requirement and soil penetration capacity of a tool, respectively (Godwin, 2007). For the models with particles of smaller nominal radii, there was no statistical differences detected in simulated soil disturbance behaviors (m, hs and P) over the entire range of PSD examined due to high standard deviations or good agreement with that measured (Zhang & Chen, 2017). In contrast, for the models with particles of larger nominal radii, all simulated soil disturbance behaviors were significantly affected by PSD (p<0.05) and had a wide range of relative errors as compared with measured values. These implied that more realistic soil loosening quality and longitudinal soil failure and accumulation degree of soil on the surface can be achieved when an appropriate PSD was employed in the DEM model. The final positive P value implied that the hardpan was more or less loosened by the subsoiler (Fig. 9). A work from Zeng et al. (2017) showed that soil porosity under the impact of a deep tillage tool initially increased and then gradually decreased to a stable value, which is in agreement with the results of the current study.

Larger soil cutting forces generally give larger soil particle forces. As shown in Fig. 13, the compressive forces of soil particles were much lower when the cutting share of the subsoiler was in the middle of the soil bin than these when the cutting share was going through the wall in the DEM simulations. The studies of Li et al. (2016) and Yang et al. (2018) indicated that soil cutting forces were much larger at the end of the travel of a soil engaging tool in the DEM simulation. The material of the physical wall of virtual soil bins was steel in this study. To avoid the abnormal rapid increase of soil cutting forces, a larger simulation bin may be practical. The higher draft force with the increase of PSD can be attributed to the larger internal friction angle of soil particles which resulted in the larger shear stress of soil particles, this is in line with the work of Coetzee & Els (2009). The predicted soil cutting forces in real time associated with the average draft and vertical forces were further examined. The characteristics of force curves of models with particles of nominal radii of 11, 13, 15 mm are similar; curves of draft and vertical forces from models with particles of nominal radii of 5, 7, 9, and 11 mm are therefore shown in Fig. 14, for the sake of conciseness. The predicted results of both draft and vertical forces were noisy and the variation was greater at larger PSDs for any nominal radii (Rn) examined. Additionally, the simulated draft and vertical forces were generally noisier at larger Rn. This explained higher standard deviations of both draft and vertical forces with the increase of nominal radii. The work of Ucgul et al. (2014b) indicated that the variation of the simulated tillage forces in real time in the DEM simulations was much greater than that measured, this agrees with the results of this study; moreover, their work also showed that the contact model between particles greatly affected both the direction and the variation of soil cutting forces in real time; in contrast, we found that soil particle-size distribution mainly affected the variation of soil cutting forces in real time. An appropriate contact model was selected by comparing the predicted and experimental vertical and draft forces in real time in the work of Ucgul et al. (2014b). This implied that more realistic DEM simulations can be gained when the variation between simulated and experimental soil cutting forces in real time is smaller.

In conclusion, for the given soil, both soil disturbance and cutting forces were significantly affected by PSD for models with particles of 9 mm nominal radii or larger; more realistic soil disturbance characteristics and draft force can be achieved by taking PSD into account in the calibration of DEM models. The significance and mechanism for the effect of PSD on soil disturbance characteristics and cutting forces provided in this study can be used to guide the decision of particle-size distribution in the soil-tool interaction DEM models (using the HMB contact model). To improve the results future work will need to consider the significance and mechanism for the effect of PSD on soil-tool interactions in different contexts (e.g. contact models between particles).

FigureFigure 12.The draft force (a) and vertical force (b) as affected by particle-size distribution (PSD) when different particle nominal radii are used in the virtual soil bin. R0: laboratory soil bin data; R5-R15: particle nominal radii ranging from 5 to 15 mm; different letters (a-e) mean significant difference at p < 0.05 for a given particle nominal radius; error bars represent standard deviations from the replicate.

FigureFigure 13.Compressive forces of soil particles: a) the moment when the cutting share of the subsoiler was in the middle of the soil bin and b) the moment when the cutting share of the subsoiler was going through the wall.

FigureFigure 14.Fluctuation of draft and vertical forces as affected by particle-size distribution (PSD) when different particle nominal radii (Rn) are used in the virtual soil bin: a) Rn = 5 mm, b) Rn = 7 mm, c) Rn = 9 mm, d) Rn = 11 mm and e) lab soil bin results.