*Review* **Review of Discrete Element Method Simulations of Soil Tillage and Furrow Opening**

**Kojo Atta Aikins 1,2,\*, Mustafa Ucgul 3,4, James B. Barr 5, Emmanuel Awuah 6, Diogenes L. Antille 2,7, Troy A. Jensen <sup>2</sup> and Jacky M. A. Desbiolles <sup>4</sup>**


**Abstract:** In agricultural machinery design and optimization, the discrete element method (DEM) has played a major role due to its ability to speed up the design and manufacturing process by reducing multiple prototyping, testing, and evaluation under experimental conditions. In the field of soil dynamics, DEM has been mainly applied in the design and optimization of soil-engaging tools, especially tillage tools and furrow openers. This numerical method is able to capture the dynamic and bulk behaviour of soils and soil–tool interactions. This review focused on the various aspects of the application of DEM in the simulation of tillage and furrow opening for tool design optimization. Different contact models, particle sizes and shapes, and calibration techniques for determining input parameters for tillage and furrow opening research have been reviewed. Discrete element method predictions of furrow profiles, disturbed soil surface profiles, soil failure, loosening, disturbance parameters, reaction forces, and the various types of soils modelled with DEM have also been highlighted. This pool of information consolidates existing working approaches used in prior studies and helps to identify knowledge gaps which, if addressed, will advance the current soil dynamics modelling capability.

**Keywords:** calibration; DEM contact models; soil dynamics; soil failure; soil forces; cohesive and frictional soils

#### **1. Introduction**

In mechanized agriculture, the energy use for soil tillage operations can be as high as 50% of the total energy used in crop production [1,2]. The energy-use efficiency associated with tillage can be increased by improving the design of tillage implements and through their correct operation and settings [3–5]. The design optimization of tillage tools and furrow openers conventionally relies on repeated prototyping and evaluation through soil bins and field experimentation. This task is laborious, time-consuming, and expensive [6,7]. In order to reduce the resource intensity involved, various analytical and numerical models for predicting soil–tool interaction and soil forces have been developed. Analytical models are typically based on the universal earthmoving equation (Equation (1)) [8,9].

$$P = (\gamma d^2 N\_\gamma + cd N\_c + c\_a d N\_a + q d N\_q) w \tag{1}$$

**Citation:** Aikins, K.A.; Ucgul, M.; Barr, J.B.; Awuah, E.; Antille, D.L.; Jensen, T.A.; Desbiolles, J.M.A. Review of Discrete Element Method Simulations of Soil Tillage and Furrow Opening. *Agriculture* **2023**, *13*, 541. https://doi.org/10.3390/ agriculture13030541

Academic Editor: Tao Cui

Received: 7 February 2023 Accepted: 14 February 2023 Published: 23 February 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

where *P* is soil cutting force (N), *γ* is the specific weight of soil (N m<sup>−</sup>3), *d* is working depth (m), *c* is soil cohesion (Pa), *ca* is soil–metal adhesion (Pa), *q* is surcharge stress (Pa), *w* is tool width (m), and *Nγ*, *Nc*, *Na*, and *Nq* are dimensionless *N* factors that are dependent on gravity, cohesion, adhesion, and surcharge, respectively.

Analytical models require a preliminary assumption of soil failure patterns to predict soil and implement forces [10,11]. These models use simple and approximate geometric profiles such as wedges and crescents that are easily expressed mathematically to model soil failure patterns [4,8,9,12,13]. They also aim to predict the maximum soil reaction force at incipient soil failure, rather than the average soil forces during the tillage process. However, in practice, soil failure occurs in more complex patterns. Analytical models also regard soil failure as bulk soil movement without accounting for interactions between individual soil particles, aggregates, and macro-organic matter evident in field conditions [7].

Numerical models, such as finite element modelling (FEM) and computational fluid dynamics (CFD), have also been employed to predict soil failure patterns and soil forces with some degree of accuracy [14,15]. However, their use is limited because they consider soil as a continuum body rather than being made up of discrete particles. Thus, continuum numerical models fail to account for variations in soil structure, physical conditions, flow, and the mixing and translocation of soil particles [6]. To overcome the shortcomings of analytical and continuum numerical methods, the discrete element method (DEM) can be employed. DEM is a discontinuum numerical method used for modelling the mechanical behaviour of granular materials. Initially developed by Cundall and Strack [16], DEM is used to model particle interactions within granular materials such as gravel, grain, soil, and powder. Unlike analytical and other numerical (continuum) approaches, DEM accounts for the discrete nature of granular particles and their interactions with neighbouring particles and interfacing objects such as contact walls and machine parts. The DEM approach involves the dynamic creation and breaking of bonds between contacting particles [7,16] and can simulate soil–tool reaction forces, as well as track particle movements and model mixing and translocating processes.

Several studies have shown that DEM can be successfully employed to model 2- and 3-dimensional space interactions between granular particles and machine parts. Operations such as grain flow in a hopper, soil movement in bulldozing, and soil deformation and displacement during field traffic, tillage and furrow opening processes have been modelled in DEM, achieving results that closely agree with experimental observations [17–21]. Figure 1 shows a comparison between a tine furrow opener interaction with soil (left) and its simulation with DEM (right) [21]. Other applications of DEM include modelling of cohesionless and cohesive/adhesive particles [10,18,20]. Discontinuum approaches such as DEM offer greater potential for more accurate prediction of soil–tool interactions in soil dynamics applications compared with continuum approaches. Thus, DEM is a powerful tool that can accurately guide and speed-up the design optimization process by researchers, developers, and manufacturers.

The objective of the work reported in this article was to review the various aspects of the application of DEM in the field of soil dynamics by focusing on soil tillage and furrow opening research for tool performance optimization. Contact models for soil particles; DEM particle size and shape considerations; calibration techniques for determining accurate input parameters; predictions of soil failure, particle movement and reaction forces; and the types of soils modelled with DEM have been reviewed. Knowledge gaps that need attention in future research, and that will go some way to advance soil dynamics modelling capability, have also been identified.

**Figure 1.** Discrete element method simulation of a narrow tine furrow opener operating in a moist sandy-loam soil. Modified from Barr et al. [21].

#### **2. Modelling Agricultural Soils with DEM**

The discrete element method is used to model soil as a collection of a finite number of individual spherical particles that interact with neighbouring particles and machine parts when subjected to external forces and forced displacement at the soil–tool interface, such as from a soil tillage tool [7]. This process induces the relative motion of particles within the bulk. Contact forces between these particles and their resultant motion are calculated using Newton's Second Law [22]. These calculations involve repeating the same algorithm at each time step of the simulation process, using results from previous calculations as input.

#### *2.1. DEM Contact Models*

DEM contact models are developed to describe the mechanical and physical interactions of granular particles with neighbouring particles or external objects. The interactions are modelled using equations of motion and contact models expressed as linear, adhesive, and elastoplastic normal contact models, as well as viscosity, tangential force, and torque models. The tangential force and torque models account for friction, rolling, and torsion [23,24]. Physical interactions between particles are expressed via combining functional elements of springs, dampers, and tangential friction. Total contact forces are expressed as the sum of spring (*F<sup>s</sup>* ) and damping (*Fd*) forces. Some commonly used contact models in DEM simulations are listed in Table 1 and reviewed below. These contact models have been implemented in commercially available software such as Bulk Flow Analyst™, Chute Analyst™, Chute Maven®, DEMpack™, Altair® EDEM™, ELFEN, GROMOS-96, ITASCA PFC (2D & 3D), LiGGGHTS®, MIMES, PASSAGE/DEM, Rocky, SimPARTIX®, StarCCM+, UDEC, 3DEC, and YADE [25,26]. Table 1 also shows the various DEM software that have been used in tillage and furrow opening research and the types of soil the contact models have been used to model.


*Agriculture* **2023** , *13*, 541

The movement of particles due to the contact forces are governed by Newton's equation of motion for linear and angular motion as expressed by Equations (2) and (3). By solving these equations, the motion of the particles can be determined. For two spherical particles (*i* = 1,2) of masses *mi* and radii *ri*, located at *xi* in contact, taking *F* as contact force, *g* as acceleration due to gravity, *Ii* as the moment of inertia of a particle, *ω<sup>i</sup>* as the angular velocity of a particle, and *Ti* as torque due to the tangential component of the contact force:

$$m\_i \ddot{\mathbf{x}}\_i = F\_i + m\_i \mathbf{g} \tag{2}$$

$$I\_i \dot{\omega}\_i = T\_i \tag{3}$$

#### 2.1.1. Linear Spring Contact Model

This contact model is linearly elastic and is the simplest contact model often used to simulate soil–tool interactions (Table 1) [6,27]. A contact force is created between the two spherical particles in contact as described above. The contact force can be decomposed into normal and tangential force components. The overlap at the contact point generates a repulsive contact force and energy dissipation. When an overlap *δ<sup>n</sup>* > 0 is formed between the two particles at a relative velocity . *δ<sup>n</sup>* in a direction normal to the contact surface, a normal contact force *Fn* is created based on the spring and dashpot models (Figure 2) such that: .

**Figure 2.** A schematic diagram of the normal (**left**) and tangential (**right**) components of the linear spring contact model.

The parameter *kn* is the normal stiffness, while *dn* is the damping coefficient. Considering an imaginary rod of radius *r* = (*r*<sup>1</sup> + *r*2)/2 between the centres of the two particles and Young's Modulus *E*:

$$F\_{\rm II} = \pi \text{Er/2} \tag{5}$$

When the tangential component of the contact force, *Ft* > *μtFn*, sliding friction occurs. The local friction coefficient *μ<sup>t</sup>* = *tan*∅, where ∅ is the internal friction angle between the particles.

The tangential component of the contact force is also given by:

$$F\_t = -k\_t \delta\_t - d\_t \dot{\delta}\_t \tag{6}$$

where *kt*, *<sup>δ</sup>t*, *dt*, and . *δ<sup>t</sup>* are tangential components of stiffness, overlap, damping coefficient, and relative velocity.

#### 2.1.2. Hertz–Mindlin Contact Model

The Hertz–Mindlin contact model (HMCM), especially when it is used with the parallel bond model (PBM, see Section 2.1.4), is the most popular contact model used by researchers [7,30,33–35,48] to simulate soil–tool interaction in tillage research. However, as an (non-linear) elastic contact model, it fails to predict vertical soil forces accurately (Table 1) [49]. In this model, the contact force consists of a non-linear Hertz component described by the hysteretic spring force-displacement relationship shown in Figure 3 and a damping component (second part of Equation (7)). It is also resolved into normal and tangential components. The HMCM and its parameters are described in Equation (7) to (16) [20,29].

**Figure 3.** Hysteretic spring force –displacement relationship used in the Hertz–Mindlin contact model (HMCM). Redrawn from Ucgul et al. [29]. The dashed line contrasts the basic linear elastic model relationship.

Normal contact force, *Fn*:

$$F\_n = -k\_n \delta\_n^{3/2} - d\_n \delta\_n^{1/4} \dot{\delta}\_n \tag{7}$$

$$k\_n = 2E\_{cq} \sqrt{r\_{cq} \delta\_n} \tag{8}$$

$$d\_n = \frac{\ln e}{\sqrt{\ln^2 e + \pi^2}} \sqrt{k\_n m\_{eq}} \tag{9}$$

where:

Equivalent radius,

$$r\_{eq} = \left(\frac{1}{r\_1} + \frac{1}{r\_2}\right)^{-1} \tag{10}$$

Equivalent Young's modulus,

$$E\_{cq} = \left(\frac{1-\nu\_1^2}{E\_1} + \frac{1-\nu\_1^2}{E\_2}\right)^{-1} \tag{11}$$

Equivalent particle mass,

$$m\_{eq} = \left(\frac{1}{m\_1} + \frac{1}{m\_2}\right)^{-1} \tag{12}$$

Tangential contact force, *Ft*:

$$F\_t = -k\_t \delta\_t - d\_t \delta\_t^{1/4} \dot{\delta}\_t \tag{13}$$

$$k\_t = 8G\_{\epsilon\eta} \sqrt{r\_{\epsilon\eta} \delta\_n} \tag{14}$$

$$d\_{l} = \frac{\ln e}{\sqrt{\ln^2 e + \pi^2}} \sqrt{k\_{l} m\_{cq}} \tag{15}$$

Equivalent shear modulus,

$$G\_{c\eta} = \left(\frac{2 - \nu\_1}{G\_1} + \frac{2 - \nu\_2}{G\_2}\right)^{-1} \tag{16}$$

#### 2.1.3. Hysteretic Spring Contract Model

The hysteretic spring contact model (HSCM) is an elastic–plastic contact model that accounts for the plastic deformation that occurs during the loading and unloading of soil. It makes the particles behave as though they undergo plastic deformation after the load reaches a yield point, as shown in Figure 4 [50]. The main disadvantage of this contact model is that it requires a large number of input parameters, making its setup and calibration of DEM material properties complex (Table 1) [49]. The HSCM comprises two parts: the spring characteristic illustrated in Figure 2 and damping. A comparative study by Ucgul et al. [29] revealed that the HSCM could model soil–tool interaction more accurately than the HMCM. The HSCM has been used to predict soil reaction forces as well as furrow profiles successfully, especially with the linear cohesion model [39–41,43]. The governing equations of the HSCM are described in Equation (17) to (20).

**Figure 4.** Hysteretic spring force-displacement relationship used in the hysteretic spring contact model (HSCM). Redrawn from Ucgul et al. [29].

Normal contact force, *Fn*:

During loading,

$$F\_n = -k\_1 \delta\_n - n\_c \sqrt{\frac{4m\_{c\eta}k\_1}{1 + \left(\frac{\pi}{\ln c}\right)^2}} \cdot \dot{\delta}\_n \tag{17}$$

During unloading and reloading,

$$F\_n = -k\_2(\delta\_n - \delta\_0) - n\_c \sqrt{\frac{4m\_{c\eta}k\_1}{1 + \left(\frac{\pi}{\ln c}\right)^2}} \cdot \dot{\delta}\_n \tag{18}$$

During unloading again,

$$F\_n = 0 - n\_c \sqrt{\frac{4m\_{c\eta}k\_1}{1 + \left(\frac{\pi}{\ln \varepsilon}\right)^2}} \cdot \dot{\delta}\_n \tag{19}$$

where *k1* and *k2* are the loading and the unloading stiffnesses, respectively, and *e* is the coefficient of restitution of the particles, and they are related as *<sup>e</sup>* <sup>=</sup> <sup>√</sup>*k*1/*k*2.

Tangential contact force, *Ft*:

$$F\_l = -n\_k k\_1 \delta\_l - \sqrt{\frac{4m\_{c\eta} n\_k k\_1}{1 + \left(\frac{\pi}{\ln c}\right)^2}} \cdot \dot{\delta}\_l\tag{20}$$

where *nk* is the stiffness factor equal to the tangential stiffness ratio to normal loading stiffness.

#### 2.1.4. Accounting for Cohesion with DEM Contact Models

In reality, agricultural soils exhibit varying levels of cohesion between particles and adhesion to walls and tools they come in contact with. Attractive pressure (that is, cohesive and adhesive forces) are induced due to the capillary effect and water bridge that exists between particles in unsaturated soils [25,51,52]. Thus, a more realistic contact model for agricultural soils should account for cohesion and adhesion. The linear cohesion and parallel bond models have been used in the DEM modelling of agricultural soils (Table 1).

#### Linear Cohesion Model

When this model is used, a cohesive or adhesive force is added to the normal force component of the contact model used for cohesionless soils. Even though the linear cohesion model itself does not include a tangential component, its addition increases the normal force, which consequently increases frictional force for greater resistance to slippage [41,50,52,53]. The linear cohesion model can be added to any of the three contact models discussed in Section 2.1.1 to Section 2.1.3 above [50]. If *Fca* (Equation (21)) is the cohesive or adhesive force, then the normal contact force is modified, as shown in Equation (22).

$$F\_{\mathfrak{C}a} = r\_{\mathfrak{C}}^2 \pi \mathfrak{C} \tag{21}$$

$$F\_n = F\_n^s + F\_n^d + F\_{ca} \tag{22}$$

$$r\_{\mathbb{C}} = \left(\frac{3r\_{\text{eq}}F\_{\text{n}}^{\text{s}}}{4E\_{\text{eq}}}\right)^{\frac{1}{3}}\tag{23}$$

The parameter *rc* is the contact radius between particles and can be determined using Equation (23). Equation (21) is called the constant cohesion model because the cohesive stress *c*ˆ is a constant. The constant cohesion model makes the model particles too sticky [52]. A modification has therefore been proposed, depending on the degree of compression

between two adjacent particles. If compression between two adjacent particles is given by Equation (24), then the cohesive stress increases with time *t* according to Equation (25).

$$
\sigma\_n = \frac{F\_n}{\pi r\_c^2} \tag{24}
$$

$$\varepsilon^{(t=n)} = k\_c \max\left(\sigma\_n^{(t=1)}, \sigma\_n^{(t=2)}, \dots, \sigma\_n^{(t=n-1)}\right) \tag{25}$$

Parallel Bond Model

An adaptation of the Hertz–Mindlin contact model (HMCM) for cohesive soils is the parallel bond model (PBM) developed by Potyondy and Cundall [54]. The PBM, based on beam theory, uses rectangular (2D) or cylindrical (3D) cement entities as parallel bonds at the point of contact between the two cohesive particles (Figure 5). This bond is modelled as an elastic beam whose length approaches zero and could be represented by a set of springs uniformly distributed over the contact plane and centred at the contact point [54]. After bond formation, normal and tangential bond forces and moment are calculated in addition to contact forces [50,55,56]. Thus, the bond can withstand or transmit both forces and moments between particles. The bond breaks when its predefined maximum normal or shear strengths are exceeded [57,58]. When no bond exists between particles, the PBM reverts to the HMCM [26]. The PBM is able to model clod formation and the brittle nature of agricultural soils in a more realistic manner [30,59]. It can be used only for particleparticle bonding, not particle-wall (tool) bonding [30,50]. The PBM is the most used model in cohesive soil tillage research [7,30,31,51,55–70]. Because the base contact model of the PBM is the HMCM, it also fails to predict vertical soil forces accurately as revealed in Table 2 [55,69]. Details of the PBM can be found in Potyondy and Cundall [54].

**Figure 5.** The parallel bond model depicted as cylindrical bond between two particles.

Johnson–Kendall–Roberts Cohesion

Another cohesion contact model combined with the HMCM is the Johnson–Kendall– Roberts (JKR) cohesion model [71]. It has been used by researchers such as Cheng et al. [33], Du et al. [36], and Zhai et al. [37] to model cohesive soils in soil tillage simulations with DEM (Table 1). Using this model, the tangential contact force and the normal and tangential damping forces of the HMCM are maintained, while the normal contact force is modified to include cohesion [36]. This modification enables the modelling of strongly adhesive bonds such as exist in dry powders or wet materials (e.g., wet soil). It captures the influence of Van der Waals forces due to contact between two surfaces [50]. A cohesion or adhesion parameter called surface energy is introduced. When this surface energy is zero, the model reverts to the HMCM. The normal contact force in HMCM-JKR is given by Equations (26) and (27), as follows:

$$F\_n = \frac{4E\_{eq}a^3}{3R\_{eq}} - 4\sqrt{\pi\gamma E\_{eq}a^3} \tag{26}$$

$$\delta\_{\pi} = \frac{a^2}{R\_{eq}} - 2\sqrt{\frac{\pi\gamma a}{E\_{eq}}} \tag{27}$$

where *a* is JKR contact radius and *γ* is surface energy (J/m2).


#### *Agriculture* **2023** , *13*, 541


**Table 2.** *Cont.* Edinburgh Elasto-Plastic Adhesion Model

Consolidation stress history is one of the main sources of cohesion in cohesive granular materials, and it must be accounted for to accurately model such materials in DEM [78]. The Edinburgh elasto-plastic adhesion model (EEPA) contact model uses a non-linear hysteretic spring model to account for the elastic–plastic contact deformation and an adhesive or cohesive force (pull-off strength) component—acting between dissimilar or similar materials, respectively based on the assumption that this force increases with increasing plastic contact area [50,79]. This model is versatile because, depending on its input parameters, it can be used as either a linear spring model (Figure 6a) or a nonlinear Hertzian spring model (Figure 6b) [79]. Figure 6 shows "a schematic diagram of particle contact and normal force-overlap (*fn*-*δ*) curve" for the EEPA contact model. A full description of the EEPA contact model can be found in Morrissey et al. [78]. The EEPA contact model has been used recently for modelling the interaction between tillage tools and agricultural soils [44–46].

**Figure 6.** Edinburgh elasto-plastic adhesion model (EEPA) normal contact force-displacement relationship (**a**) linear and (**b**) non-linear [50].

#### *2.2. Particle Size and Shape*

Particle shape and size used in DEM significantly affect the necessary simulation time and the accuracy of simulation results [10,28,38]. They are input parameters that should be carefully chosen during calibration for DEM particle interactions to be as close to reality as possible [20,49]. In DEM simulations, it is ideal to use particles of similar sizes to the actual granular materials being modelled. For instance, the actual sizes of agricultural soil particles are relatively small, ranging from several nanometres to about 2 mm for very coarse sand [80]. To model actual particle sizes in DEM requires unrealistically long computation time and impractically high computer processing power [78]. The most timeconsuming part of soil particle DEM simulations is contact detection, and is proportional to the number of particles [81]. For this reason, larger particle sizes than real soil particle sizes are generally employed in DEM [20,38,76]. The larger particles are sometimes implied to represent soil aggregates instead of individual soil particles and somewhat capture the bulk behaviour of a structured soil profile [78].

In reality, soil particles come in various irregular shapes. Thus, particles used in DEM should be not only of a similar size range but also of a similar shape range to actual soil particles to ensure simulations are more representative of realistic bulk behaviour. The basic shape of a DEM particle is a sphere (or circle in 2D modelling) under most DEM codes [20]. Spherical particles approximate and simplify simulations by improving contact detection efficiency and reducing computation time [82]. Usually, irregular (non-spherical) particles are created by clumping a number of spherical particles together, as shown in Figure 7 [28,32,46,83]. This enables the use of spherical particle contact detection algorithms that are simpler and require shorter computation time than those of irregular shapes [24,84].

**Figure 7.** Different associations of spherical DEM particles used to more realistically account for the effects of soil particle shape on bulk behaviour.

Nonetheless, clumped particles also require relatively higher computational time than when purely spherical particles are used. Therefore, most studies adopt spherical particles to represent soil particles or soil aggregates in DEM simulations. Spherical particle assemblies simulating the soil profile are characterized by considerably lower internal friction and shear strength than actual soil particles due to lower impact of rolling friction. However, this is usually overcome by introducing an arbitrary high rolling friction coefficient to simulate the interlocking tendencies that exist between the irregular shape soil particles [83,85].

#### **3. Calibration Techniques for Determining DEM Input Parameters**

Running a DEM model involves providing it with such input parameters aimed at simulating soil behaviour as close to real soils as possible. Accurate results can only be obtained with accurate input parameters [57,86]. Several approaches exist for calibrating DEM input parameters that accurately represent both soil to soil particle properties and soil particle to tool or machine interface properties. The most common calibration methods for the former include the angle of repose and hopper discharge, direct shear and triaxial tests, and corresponding in situ soil measurements. The most common calibration methods for the latter include the inclined plane test, the modified shear test, and corresponding in situ measurements. All these approaches are focused on bulk responses (i.e., natural stable state and force reactions) of soil under an applied load. After experimental runs in the laboratory or field, these experiments are then replicated numerically as closely as possible, optimizing parameters iteratively until bulk numerical responses agree with field or laboratory measurements [20]. Trial and error methods have traditionally been relied upon in the past while, more recently, the application of response surface methodology (RSM) is demonstrating benefits of significantly reducing the number of numerical simulations required for accurate calibration [34,45,87–91].

#### *3.1. Angle of Repose Test*

The angle of repose test is used to assess flowability and inter-particle friction of loose soil [92–94]. This test is also essential when there is a need for a qualitative assessment of soil surface and furrow profiles in tillage simulations [29]. Various researchers [21,42,84,88,95] used the angle of repose test to calibrate coefficients of static and rolling friction between soil particles.

In this test, the soil is allowed to flow by gravity onto a flat surface to form a cone pile. The angle of repose is measured as shown in Figure 8a [29]. Another approach for the angle of repose test is to confine the particles being modelled within the walls of a box, ensuring the top of the particles is levelled. By removing one of the sidewalls, the particles flow to form the angle of repose as shown in Figure 8b [83]. This test is usually used for cohesionless particles and particles with low cohesion with good flowability. However, the general principle is that repeatable observations can be made during key stages of the "angle of repose" experiment. For instance, Roessler and Katterfeld [96] reported successfully calibrating DEM parameters for cohesive soil using the angle of repose test. A cylinder was filled with the cohesive soil, and the cylinder was gradually lifted as shown in Figure 8c. Reproducible phases of soil flow were observed, namely "the build-up of a stable bulk material column, the convex bending of the column, and the beginning of collapse of the column." Aikins et al. [41] observed a reproducible dome-like pile of cohesive soil (Figure 9a,b) and used the results to calibrate soil–soil coefficients of static and rolling friction values.

**Figure 8.** Angle of repose test setup used by (**a**) Ucgul et al. [29] for a cohesionless soil, (**b**) Mousaviraad et al. [83] for maize, and (**c**) Roessler and Katterfeld [96] for a cohesive soil.

**Figure 9.** Repeatable dome-like soil pile obtained in the angle of repose test for a loose cohesive soil by Aikins et al. [41]: (**a**) laboratory experiment and (**b**) optimized DEM simulation.

#### *3.2. Inclined Plane Test*

The inclined plane test has been used to determine soil–tool or soil–machine coefficients of static and rolling friction [29]. A schematic diagram of the setup for the inclined plane test is shown in Figure 10. A flat bed of the soil to be modelled is packed into a tray and held on a table with adjustable horizontal inclination. A block of tool material and ball bearings are separately placed on the flat bed, and the table is tilted to an angle Ψ at which the block just starts sliding or the ball just starts rolling down the inclination. The block is used for the determination of the soil–tool coefficient of static friction (*μs*), while the ball bearing is used in the determination of the soil–tool coefficient of rolling friction (*μr*). If the mass of the block is *ms*, the mass of the ball is *mr*, and the angles at which sliding and rolling occur are Ψ*<sup>s</sup>* and Ψ*r*, respectively, then the coefficients are calculated according to Equations (28) and (29) [29,97].

**Figure 10.** Schematic diagram of inclined plane test setup with a block [97] for static friction coefficient. A ball can be used as well for rolling friction coefficient as found in Ucgul et al. [29].

Soil–tool coefficient of static friction,

$$\mu\_s = \frac{m\_s g \sin \Psi\_s}{m\_s g \cos \Psi\_s} = \tan \Psi\_s \tag{28}$$

Soil–tool coefficient of rolling friction,

$$\mu\_{\varGamma} = \frac{m\_r g \sin \Psi\_r}{m\_r g \cos \Psi\_r} = \tan \Psi\_r \tag{29}$$

#### *3.3. Direct Shear Test*

The direct shear test is used to determine internal soil parameters namely, cohesion and internal friction angle (for soil-to-soil particle interactions). The modified shear test is used to determine the adhesion and external friction angle (for soil to tool or machine interface properties). These are typically used as direct DEM input parameters to support the calibration of other parameters. It has been used for cohesive and adhesive soils, as well as cohesionless soils [41,44,53].

This approach uses the normal and shear stresses acting on a column of granular materials' cross-section. The experimental setup consists of two shear boxes, one placed on the other and filled with the granular material being modelled. One half is fixed while the other is made movable horizontally in one direction (Figure 11a). A specified normal force (*Fa*) is applied while an increasing horizontal (shearing) force (*Fb*) is applied to the movable half till a certain amount of displacement occurs [98]. At that point, the horizontal force would have reached a maximum value and remain constant or slightly increase or decrease afterward [30]. The experiment is repeated several times with different normal force values.

**Figure 11.** (**a**) Direct shear box schematic, and (**b**) normal and shear stress expected relationship.

In the modified shear box test, the bottom half of the shear box is replaced with a material matching that of the tool or machine part. Normal force and shearing forces are applied the same way as described above. Corresponding normal and shear stresses are plotted as shown in Figure 11b. Given that the cohesive strength of the soil is *c*, the cross-sectional area of the shear box is *A*, and the internal friction angle is *φ*, the relationship between the normal force and the horizontal force is displayed in Equation (30).

$$F\_b = cA + F\_a \tan \phi \tag{30}$$

#### *3.4. Triaxial Compression Test*

This test can also determine the soil cohesive strength and internal friction angle and used to calibrate DEM input parameters for cohesive-frictional soils [52,72]. A triaxial compression test consists of loading an undisturbed cylindrical soil specimen insulated with an impermeable membrane and subjected to an adjustable confining pressure within a water chamber [18,99,100]. The specimen is then subjected to combined axial (σ1) and radial (σ3) stresses as indicated in Figure 12 until soil failure is achieved [98]. The radial stress (σ3) is first applied around the specimen to a set level via the confining water pressure. An axial strain is then mechanically applied at a controlled rate which generates a corresponding additional deviator stress (*q*) logged over time and combining with σ<sup>3</sup> to form a resultant axial stress σ1. The above steps are repeated several times under increasing radial stress. The plots of the deviator stress (*q* = σ<sup>1</sup> – σ3) against axial strain identify each deviator stress value at failure and a simple process—for instance using the Mohr circle method—is then used to quantify soil cohesion and internal friction angle [101].

**Figure 12.** Specimen stress state in a triaxial compression test: retrieved from (**a**) Rees [99] and (**b**) Ahlinhan et al. [102].

#### *3.5. In Situ Approaches*

Measurements of soil mechanical properties are most accurately done in the laboratory [103]. However, while laboratory methods can precisely measure soil properties, samples may not always be fully representative of field soil conditions due to sampling and handling limitations and time-related changes between sampling and testing. Hence, some researchers have used in situ approaches to measure soil properties for DEM calibration purposes. Kim et al. [44] used an on-site measurement system to determine soil mechanical properties such as shear modulus, Young's modulus, and soil–tool static and rolling friction (Figure 13). On the other hand, Aikins et al. [41] used an on-site cone penetration test (Figure 14) to calibrate Young's modulus of a Black Vertosol.

**Figure 13.** Instruments for in situ measurements of soil mechanical properties: (**a**) sampled undisturbed soil specimen, (**b**) soil–metal static friction, (**c**) soil–metal rolling friction, (**d**) direct shear test, (**e**) shear modulus, and (**f**) Young's modulus. Retrieved from Kim et al. [44].

**Figure 14.** (**a**) Motorised cone penetrometer for in situ measurements, and (**b**) DEM cone penetration simulation used by Aikins et al. [41].

Asaf et al. [10] proposed grouser shear and sinkage or penetration tests using wedges of different wedge angles and a plate for calibrating DEM contact parameters. Jang et al. [82] also used a rectangular plate, while Ucgul et al. [38] and Ucgul et al. [29] used circular disc and cone penetration tests to calibrate model parameters. Cheng et al. [33] used a soil adhesion mass test to determine DEM input parameters of wet clay soil by employing the Plackett–Burman test and response surface methodology (RSM) to optimise input parameters.

#### **4. Prediction of Soil Failure, Loosening, and Disturbance Parameters**

*4.1. Soil Failure and Loosening*

Tamas et al. [30] and Barr et al. [21] have revealed the ability of DEM to predict soil rupture and crack propagation, which is an advantage over FEM. Some researchers have used velocity profiles [7,32] or displacement profiles [41,69,104] as soil loosening indicators. Others [21,30] used porosity (in PFC3D Particle Flow Code) or voidage (in EDEM 2.7TM), respectively, to measure the degree of particles loosening in DEM. In the work of Tamas et al. [30], for instance, it was found that the DEM modelled soil porosity and soil-break-up resulting from loosening by sweeps increased with both increasing speed and rake angle, which agrees with experimental results.

Identifying particle movement or loosening is mainly used in defining the boundary between disturbed and undisturbed particles to simulate soil failure boundary or furrow profile. Barr et al. [21] argued that using velocity and displacement profiles is based on the assumption that particle movement results in only soil loosening, ignoring the fact that particle movement also occurs during a soil compaction process. The validity of this assumption is therefore limited to tools operating above their critical depth. Additionally, these approaches are open to subjective decisions since a threshold has to be arbitrarily defined to differentiate between the "so-called" loosened and unloosened particles. For example, Murray [69] had to describe loosened particles as having a displacement magnitude above 5 mm. Barr et al. [21] instead proposed and used a voidage grid (Figure 15) to define failure boundaries. A voidage grid was applied to the DEM particles after the tillage process was completed and the particles had settled, which reflects experimental practice. Voidage is similar to soil porosity as it measures the proportion of volume not occupied by particles. An empty space will have a voidage of 100%, while a completely filled space will have a voidage of 0%. With *Vg* being grid volume and *Vp* the total volume of particles whose centroids are located within the grid, voidage can be calculated according to Equation (31).

$$Voolage = \frac{V\_{\ $} - V\_{p}}{V\_{\$ }} \times 100\% \tag{31}$$

**Figure 15.** Furrow profiles created using voidage grid bins in EDEM 2.7TM. Adapted from Barr [49].

Aikins et al. [41] used particle displacement (PD) analysis to determine the loosened furrow boundary in DEM. The displacement threshold was not set arbitrarily as was done by Murray [69]. Aikins et al. [41] defined the loosened furrow boundary based on two criteria:

1. Minimum particle displacement caused directly by an opener occurs with particles just adjacent to the bottom part of the opener (for wide tines) or particles aligning the walls of the slot below critical depth (for narrow tines).

2. To establish a sharp contrast between displaced and undisturbed particles, particle locations immediately after particle loosening (i.e., before the particle settle) has to be used.

Aikins et al. [41] traced the minimum particle displacement up the profile to produce a loosened furrow boundary as shown in Figure 16a. Figure 16a is a contour plot of the width and depth of the virtual soil bin profile against displacements (resultant) for each particle within the profile. In some studies [21,32], disturbed soil surface profile after tillage was determined using voidage grid binning and velocity profile. However, Aikins et al. [41] used the profile of the top surface of displaced DEM particles as the disturbed surface profile after the particles had settled because that gives more realistic results and is similar to what actually happens in field experiments. Wang et al. [74] employed another approach using the "clipping" module in EDEM simulation software to define the disturbed soil boundary (Figure 16b). The furrow profile was obtained by connecting the boundaries of the different layers of disturbed soil.

**Figure 16.** Definitions of furrow boundary used by (**a**) Aikins et al. [41] and (**b**) Wang et al. [74].

#### *4.2. Soil Movement and Disturbance Parameters*

From the loosened criteria described in the previous section, furrow profile, soil movement, and various soil disturbance parameters have been predicted or determined in DEM simulations with varying levels of relative error (RE) compared to experimental results. Such soil disturbance parameters include the lateral, forward, and upward movement of particles; furrow width at soil surface; loosened furrow cross-sectional area; furrow % backfill and dip area. Barr et al. [21] found an RE of 9% in loosened furrow cross-sectional area, 26% in furrow width, 14% in dip area, 0.8% in furrow backfill, 16% in ridge height, and 9% in lateral soil throw in a DEM prediction of soil disturbance parameters with narrow point openers operating in a sandy-loam soil. Barr and Fielke [105] closely predicted lateral soil throw and soil layer mixing using narrow tine openers with 35◦ and 90◦ rake angles and a bentleg opener (Figure 17). Using furrow openers with different rake angles and cutting edge cross-sections, Aikins et al. [41] closely predicted furrow profiles and similar patterns for surface profiles. The majority of DEM predictions of furrow cross-sectional area, furrow width, critical depth, and lateral soil throw had an RE of 1% to 20%. However, poor predictions were made for ridge height due to the use of large DEM particles (radius of 5 mm) [41].

**Figure 17.** Comparison of furrow profiles for 16 mm wide narrow openers at different rake angles (35◦, 53◦, 72◦, and 90◦) obtained from soil bin experiment and DEM simulations. Retrieved and modified from Barr [49].

Wang et al. [74] determined soil looseness, furrow width, soil disturbance coefficient, and soil disturbance area ratio with an RE from 3.24% to 41.64% for a winged subsoiler and 0.24% to 28.74% for a non-winged subsoiler. There was also a satisfactory agreement in "shape and magnitude" of lateral, forward, and upward displacement of different soil layers resulting from a sweep cultivator between experimental and DEM simulation results [62]. Using a disc with tilt angles from 0◦ to 20◦, Murray [69] estimated an average absolute RE of about 10.53% for lateral soil throw. The DEM simulation also revealed, in agreement with the experimental result, that lateral soil throw increased with increasing tilt angle. For a hoe furrow opener, a RE of 14.8% was recorded.

Reduction in the forward movement of soil particles at greater depth has been predicted through DEM simulations [21,32,41]. Other researchers have previously described this phenomenon [9,106–108] through experimental work and analytical models in the concept of critical depth. Below the critical depth, soil movement changes from forward, sideways, and upward directions (generating a loosened crescent-shaped soil failure) to mainly forward and sideways, generating a compaction type failure in a horizontal plane. Barr et al. [21] and Aikins et al. [41] closely predicted the narrowing of furrow down the profile and critical depth with a 90◦ rake angle (Figures 15 and 16a). Hang et al. [32] observed a reduction of the forward movement of particles in the inter-row zone between tines with increasing tine spacing in both DEM and experimental results. Greater inter-row zone soil movement with narrower tine spacing was attributed to more intensive interaction between tines. DEM can also simulate greater soil movement with wing attachment to tine cutting tools [31].

Greater soil upheaval observed with a low rake angle opener (35◦) and increased lateral throw of soil (due to splashing effect) typically found with steeper (90◦) rake angle openers have been successfully replicated with DEM [21]. However, Ucgul et al. [38] observed that the dynamic height of soil flow under sweep openers was under-predicted by 23% to 35% at speeds of 5 to 12.5 km/h and did not follow the observed shape using spherical particles of 10 mm radii. The prediction was improved by using smaller particles of radii 1.5 mm, which still underpredicted soil flow height by 15%. Using particles of smaller radii, however, considerably increased computation time. Lateral soil throw was also under predicted by about 32% and 9% with radii 10 and 1.5 mm, respectively.

Chen et al. [7] estimated a maximum of 3% RE in furrow cross-sectional area, up to 4% RE in disturbed width, from 14% to 26% more soil (by volume) heaped above the soil surface and 5% to 15% more emptied cross-section below the soil surface. Overall, a close agreement was observed between experimental and simulation results. Hang et al. [32] estimated less than 20% RE between DEM predicted and experimentally determined soil disturbance coefficient and soil looseness. Saunders et al. [76] reported significant underprediction and correlations between measured and predicted furrow area (r = 0.82) and maximum soil throw (r = 0.88) when optimizing the performance of a mouldboard skimmer in a sandy-loam soil. Several furrow parameters from bentleg openers operating in sandy-loam soil, including loosened cross-sectional area (RE = 14.9%), furrow dip area (RE = 14.4%), backfill (RE = 1.8%), ridge height (RE = 16.8%), and lateral soil throw (RE = 14.9%), were accurately predicted using the voidage grid bin approach [39]. Furthermore, Barr et al. [39] observed the same findings as those measured in soil bin investigations, with DEM simulations of bentleg openers also achieving 100% backfill and cancelling furrow spill over.

The DEM has also been used to simulate rotary tiller operations for design optimization. Zhang et al. [109] and Hirasawa et al. [110] closely predicted the height and pattern of soil surface undulations after rotary tillage in DEM. Soil movement pattern during rotary tillage and improving soil layer mixing with tillage depth and travel speed for a rotary tiller have also been closely predicted [36,109]. Cheng et al. [33] recorded an RE of 1.84% in mass of soil that adhered to rotary tiller blades.

#### **5. Prediction of Tillage Forces**

Most attention in the DEM simulation of soil cutting tools has been on predicting soil forces. Table 2 shows relative error (RE) values in draught and vertical force predictions with DEM, travel speed and operating depths for the stated tillage tools, soil types, and DEM contact models. Draught force cycles between peaks at incipient soil failure and troughs at the start of the reloading phase have been well captured in DEM simulations of sweeps by Tamas et al. [30]. In the same DEM study, draught was found to increase with greater speed and rake angle. Draught was predicted with 4% to 12% RE as speed was increased from 0.5 to 2.4 m s−<sup>1</sup> with a sweep tine. Bo et al. [31] observed a similar trend between draught force measured for four subsoilers in the soil bin and that obtained through DEM simulation. A winged subsoiler among the four had the highest draught force in both soil bin tests (up to 50% more) and simulations (up to 55% more). With all four subsoilers, DEM predicted draught force with relative errors below 4%.

Additionally, Wang et al. [74] showed that a winged subsoiler operating at a speed and depth of 3 km h−<sup>1</sup> and 300 mm, respectively, had an RE of 9.71%, whereas the nonwinged subsoiler obtained an RE of 15.08% when the draught force was compared with the experimental result. Chen et al. [7] observed about 4–31% RE between the draught of experimental and DEM results. A good correlation was obtained between the measured and predicted draught forces (r = 0.95), whereas a more limited correlation was observed for vertical force (r = 0.71). With blunt (R90B) and chamfered (C2S) narrow openers with a 90◦ rake angle, a blunt opener with a 45◦ rake angle, and a bentleg opener, Aikins et al. [41] predicted draught force with REs of 20%, 22%, 31%, and 5%, respectively. Vertical force was also predicted with 8% and 20% relative error for R90B and C2S, respectively, but poorly for the two other openers.

DEM simulations with a mouldboard plough were also able to simulate the gradual entry of the mouldboard into the soil with a gradual draught increase. For the two soil conditions used in this study, 9% and 2.4% errors in cultivator tool draught were observed for a soft-wet soil and a hard-dry soil, respectively. DEM also closely predicted an increase in draught with increasing depth, with RE ranging from about 3% to 15% [56]. Kim et al. [44] observed that draught force increased with increasing depth with an overall average RE of 7.45%. Tong et al. [73] showed that, across four tillage depths, simulated draught and vertical forces were up to 10% smaller than those measured in the field.

In a DEM prediction of horizontal and vertical soil forces with DEM particles clumped to form different shapes (similar to that shown in Figure 7), Ono et al. [28] obtained the most accurate predictions with the three linearly overlapping spheres. The worst prediction was obtained with simple spherical particles. Ucgul et al. [38] observed a linear increase in draught force against sweep tine width measured experimentally and predicted using DEM with a maximum RE of 8%. Likewise, a non-linear increase in vertical force against width with a maximum RE of 13.7% was recorded. High correlations were recorded between measured and predicted draught forces (r = 0.978) and vertical forces (r = 0.971) with tool speeds from 5 to 12.5 km h−<sup>1</sup> and a depth of 70 mm. Prediction of the effect of rake angle on soil forces followed a similar trend (r values of 0.98 and 0.97) and had an RE of 11.6% and 15.2% for draught and vertical forces, respectively. Ucgul et al. [111] again obtained an accurate prediction of draught and vertical forces of a sweep tillage tool at varying speeds and geometry with r values ranging from 0.84 to 0.92. Murray [69] estimated an average RE of 1.86% for draught and 50.7% for vertical force with a flat single disc opener. For rotary tillers, Zhang et al. [109] reported a 12% RE in power consumption, while Du et al. [36] predicted increasing torque with tillage depth (150 to 180 mm) and travel speed (about 2 to 3 km h<sup>−</sup>1).

#### **6. Soils Modelled in DEM Simulations**

Tables 1 and 2 list soil types used in various tillage and furrow opener DEM simulations and their bulk densities, soil water contents, and cohesive strengths. It can be seen that most of the soils modelled with DEM are of sandy to sandy-loam textures. DEM modelling of highly cohesive soils is still relatively scarce in the literature. Bravo et al. [18] used DEM to model highly cohesive clay soil (Vertosol) with cohesive strength of up to about 125 kPa when the soil was highly compacted (bulk density of about 1400 kg m<sup>−</sup>3) and relatively dry (soil water content of about 18%). Aikins et al. [41] also modelled a Vertosol with cohesive strength of 46.4 kPa.

The properties and flow characteristics of sandy soils differ from those of clay soils, which show cohesive and adhesive properties in the presence of sufficient moisture. The successful modelling of a Vertosol and its interaction with tools in DEM by Bravo et al. [18] and Aikins et al. [41] revealed the ability of DEM to model cohesive soils. Although some authors [33,41,43,44,104,112,113] have recently modelled cohesive soils using DEM, more attention is needed in future research to cover the wide spectrum of agricultural soils.

#### **7. Conclusions**

Based on this review, the following conclusions can be drawn:


Based on the review conducted, the following recommendations are made for future research:


**Author Contributions:** Writing—original draft preparation, K.A.A.; writing—review and editing, K.A.A., M.U., J.B.B., E.A., D.L.A., T.A.J. and J.M.A.D.; visualization, K.A.A. and E.A.; supervision, M.U., J.B.B., D.L.A., T.A.J. and J.M.A.D.; project administration, D.L.A., T.A.J. and J.M.A.D.; funding acquisition, K.A.A. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by CLAAS Stiftung (Harsewinkel, Germany) and the University of Southern Queensland (Toowoomba, QLD, Australia).

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors are grateful to the Centre for Agricultural Engineering at the University of Southern Queensland (Toowoomba, Qld, Australia), the Agricultural Machinery Research and Design Centre at the University of South Australia (Adelaide, SA, Australia), CSIRO Agriculture and Food (Canberra, ACT, Australia), and the CLAAS Foundation (Harsewinkel, Germany, http://www.claas-stiftung.com/ accessed on 17 February 2023) for financial and operational support to conduct this review.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**


#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Agriculture* Editorial Office E-mail: agriculture@mdpi.com www.mdpi.com/journal/agriculture

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel: +41 61 683 77 34

www.mdpi.com ISBN 978-3-0365-7295-6