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

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.

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,4,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 = (γd²N_γ + cdN_c + c_adN_a + qdN_q)w

(1)

where P is soil cutting force (N), γ is the specific weight of soil (N m⁻³), d is working depth (m), c is soil cohesion (Pa), c_a is soil–metal adhesion (Pa), q is surcharge stress (Pa), w is tool width (m), and N_γ, N_c, N_a, and N_q 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,18,19,20,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.

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^s) and damping (F^d) forces. Some commonly used contact models in DEM simulations are listed in

Table 1. Discrete element method contact models, their advantages and disadvantages, types of soil modelled, and software used.

Contact Model	Advantages	Disadvantages	Types of Soil Modelled	References	Software Used by Researchers
Linear spring contact model	Simple to use.	Does not account for nonlinearity in loading and unloading cycles and plastic deformation of soil.	Sandy	Tanaka et al. [27], Asaf et al. [10], Shmulevich et al. [6], Ono et al. [28]	PFC2D, EDEM
Linear spring contact model with cohesion	Allows users to consider cohesion in the linear spring contact model.	Only considers the cohesion through the normal direction.	Vertosol	Bravo et al. [18]	DEMeter++
Hertz–Mindlin contact model	Simple to use. Although this model was designed for fine, dry particles, it can be used to model wet particles as well.	Inaccurate prediction for vertical tillage force.	Sandy	Ucgul et al. [29]	EDEM
Parallel bond model (PBM) or Hertz–Mindlin contact model with cohesion	Allows users to model cohesion.	Excessive forces cause the bonds to be broken irrationally.	Coarse sand, loamy, sandy loam, loessal, clay, sandy clay loam, loamy clay	Tamas et al. [30], Chen et al. [7], Bo et al. [31], Hang et al. [32], Cheng et al. [33], Yang et al. [34], Hoseinian et al. [35]	EDEM PFC3D
Hertz–Mindlin contact model with Johnson–Kendall–Roberts (JKR)	Enables the modelling of strongly adhesive bonds such as exist in dry powders or wet materials.	-	Clay, silty clay loam	Cheng et al. [33], Du et al. [36], Zhai et al. [37]	EDEM
Hysteric spring contact model	Accounts for plastic deformation during loading and unloading of soil. Suitable for both cohesive cohesionless soils.	It requires a large number of input parameters, making its setup and calibration complex.	Sandy	Ucgul et al. [38], Ucgul et al. [29]	EDEM
Hysteric spring contact model with linear cohesion contact model	Allow users to consider cohesion in the hysteretic spring contact model.	Only considers the cohesion through the normal direction.	Sandy loam, clay (Vertosol)	Barr et al. [21], Barr et al. [39], Makange et al. [40], Aikins et al. [41], Awuah et al. [42], Wang et al. [43]	EDEM
Edinburgh elasto-plastic adhesion model	It is versatile since it can be used as a linear or non-linear Hertzian spring model. It also allows tensile forces to develop and a non-linear force-displacement behaviour in compression.	It requires a lot of input parameters.	Clay, clay loam, sandy loam, loam, sandy	Kim et al. [44], Wu et al. [45], Zhao et al. [46], Sun et al. [47]	EDEM, PFC3D, LiGGGHTS

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.

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 m_i and radii r_i, located at x_i in contact, taking F as contact force, g as acceleration due to gravity, I_i as the moment of inertia of a particle, ω_i as the angular velocity of a particle, and T_i as torque due to the tangential component of the contact force:

$$m_i{\ddot{x}}_i=F_i+m_{i}g$$

(2)

$$I_i{\dot{\omega }}_i=T_i$$

(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 ${\delta }_n>0$ is formed between the two particles at a relative velocity ${\dot{\delta }}_n$ in a direction normal to the contact surface, a normal contact force $F_n$ is created based on the spring and dashpot models (Figure 2) such that:

$$F_n=k_n{\delta }_n+d_n{\dot{\delta }}_n$$

(4)

The parameter k_n is the normal stiffness, while d_n is the damping coefficient. Considering an imaginary rod of radius $r=\left(r_1+r_2\right)/2$ between the centres of the two particles and Young’s Modulus E:

$$F_n=\pi Er/2$$

(5)

When the tangential component of the contact force, $F_t>{\mu }_t{F}_n$, sliding friction occurs. The local friction coefficient ${\mu }_t=tan\varnothing$, where $\varnothing$ 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$$

(6)

where k_t, ${\delta }_t$, d_t, and ${\dot{\delta }}_t$ 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,34,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].

Normal contact force, F_n:

$$F_n=-k_n{\delta }_n^{3/2}-d_n{\delta }_n^{1/4}{\dot{\delta }}_n$$

(7)

$$k_n=2E_{eq}\sqrt{r_{eq}{\delta }_n}$$

(8)

$$d_n=\frac{\ln e}{\sqrt{{\ln }^{2}e+{\pi }^2}}\sqrt{k_n{m}_{eq}}$$

(9)

where:

Equivalent radius,

$$r_{eq}={\left(\frac{1}{r_1}+\frac{1}{r_2}\right)}^{-1}$$

(10)

Equivalent Young’s modulus,

$$E_{eq}={\left(\frac{1-{\nu }_1^2}{E_1}+\frac{1-{\nu }_1^2}{E_2}\right)}^{-1}$$

(11)

Equivalent particle mass,

$$m_{eq}={\left(\frac{1}{m_1}+\frac{1}{m_2}\right)}^{-1}$$

(12)

Tangential contact force, F_t:

$$F_t=-k_t{\delta }_t-d_t{\delta }_t^{1/4}{\dot{\delta }}_t$$

(13)

$$k_t=8G_{eq}\sqrt{r_{eq}{\delta }_n}$$

(14)

$$d_t=\frac{\ln e}{\sqrt{{\ln }^{2}e+{\pi }^2}}\sqrt{k_t{m}_{eq}}$$

(15)

Equivalent shear modulus,

$$G_{eq}={\left(\frac{2-{\nu }_1}{G_1}+\frac{2-{\nu }_2}{G_2}\right)}^{-1}$$

(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,40,41,43]. The governing equations of the HSCM are described in Equation (17) to (20).

Normal contact force, F_n:

During loading,

$$F_n=-k_1{\delta }_n-n_c\sqrt{\frac{4m_{eq}{k}_1}{1+{\left(\frac{\pi }{\ln e}\right)}^2}}·{\dot{\delta }}_n$$

(17)

During unloading and reloading,

$$F_n=-k_2\left({\delta }_n-{\delta }_0\right)-n_c\sqrt{\frac{4m_{eq}{k}_1}{1+{\left(\frac{\pi }{\ln e}\right)}^2}}·{\dot{\delta }}_n$$

(18)

During unloading again,

$$F_n=0-n_c\sqrt{\frac{4m_{eq}{k}_1}{1+{\left(\frac{\pi }{\ln e}\right)}^2}}·{\dot{\delta }}_n$$

(19)

where k₁ and k₂ are the loading and the unloading stiffnesses, respectively, and e is the coefficient of restitution of the particles, and they are related as $e=\sqrt{k_1/k_2}$.

Tangential contact force, F_t:

$$F_t=-n_k{k}_1{\delta }_t-\sqrt{\frac{4m_{eq}{n}_k{k}_1}{1+{\left(\frac{\pi }{\ln e}\right)}^2}}·{\dot{\delta }}_t$$

(20)

where n_k 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 F_ca (Equation (21)) is the cohesive or adhesive force, then the normal contact force is modified, as shown in Equation (22).

$$F_{ca}=r_c^2\pi \widehat{c}$$

(21)

$$F_n=F_n^s+F_n^d+F_{ca}$$

(22)

$$r_c={\left(\frac{3r_{eq}{F}_n^s}{4E_{eq}}\right)}^{\frac{1}{3}}$$

(23)

The parameter r_c 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 $\widehat{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}$$

(24)

$${\widehat{c}}^{\left(t=n\right)}=k_c\underset{}{\max }\left({\sigma }_n^{\left(t=1\right)}, {\sigma }_n^{\left(t=2\right)},\dots ..{\sigma }_n^{\left(t=n-1\right)}\right)$$

(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 particle-particle 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,56,57,58,59,60,61,62,63,64,65,66,67,68,69,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].

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}$$

(26)

$${\delta }_n=\frac{a^2}{R_{eq}}-2\sqrt{\frac{\pi \gamma a}{E_{eq}}}$$

(27)

where a is JKR contact radius and γ is surface energy (J/m²).

Table 2. Relative errors in DEM-predicted soil–tool reaction forces, travel speeds, and operating depths reported in the various literature across soil types.

Reference	Relative Error in DEM Prediction (%) Relative to Measured Data		Travel Speed (km h−1)	Operating Depth (mm)	Tillage Tools	Soil Texture	Dry Bulk Density (kg m−3)	Soil Water Content (%, w/w)	Cohesive Strength (kPa)	Contact Model
	Draught	Vertical Force
Sadek et al. [58]	n/a	n/a	n/a	n/a	n/a	Sandy soil	990 1280 1360 1500	0.02 13 21.5	1.23-32.70	PBM
Chen et al. [7]	4 to 31	n/a	3.19 (average)	100	Sweep tine	Coarse sand Loamy sand Sandy loam	1410 1330 1410	8.98 14.84 18.2	15.7 25.2 36	PBM
Obermayr et al. [52]	n/a	n/a	2.16–4.5	10–200	Bulldozer blade	n/a	1900	n/a	11.16	LSCM + cohesion
Tamas et al. [30]	4 to 12	n/a	1.8–8.64	200	Sweep tine	Sandy soil	1850	6.33	11.86	PBM
Bravo et al. [18]	9, 24	n/a	-	150–500	Para-plough and mouldboard plough	Clay (Vertosol)	1000 1200 1400	8 18 20 35	25–125	LSCM + cohesion
Li et al. [56]	3 to 15	n/a	3.6	180–260	Subsoiler	n/a	n/a	19	n/a	PBM
Mak and Chen [61]	n/a	n/a	2.2–6.59	50–200	Sweep tine	Loamy sand	1320	11.3	13.9	PBM
Obermayr et al. [72]	n/a	n/a		100–200	Straight-vertical blade and bulldozer blade	Sand	1520 1980 1870	10 15	6–22.5	LSCM + cohesion
Ucgul et al. [38]	≤11.6	≤15.2	5–12.5	70	Sweep tine	Sandy loam	1750	8	6	HSCM
Ucgul et al. [53]	n/a	n/a	4–12	75	Sweep tine	Sandy loam	1320 1780 1880	1 15 13	3 15 22	HSCM + LCM
Kotrocz et al. [60]	n/a	n/a	n/a	50–150	Cone penetrometer	Loamy sand	1632	15.8	6.61–8.66	PBM
Li et al. [70]			2.99	3–18	Claw	Sandy loam	1300	n/a	17.5	PBM
Murray [69]	1.86	50.7	8	38	Disc and hoe openers	Clayey lacustrine	1560	19.6	n/a	PBM
Hang et al. [32]	n/a	n/a	3	300	Subsoiler	Loamy clay	1346	12.5	11.8	PBM
Milkevych et al. [62]	n/a	n/a	3.2	100	Sweep tine	Coarse sand Loamy sand	1410 1330	9 14.8	15.8 25.1	PBM
Tekeste et al. [55]	9, 12	-59, -49	0.79–9.65	102	Sweep tine	Loam	1307	8.99	33	PBM
Tong et al. [73]	<10	<10	7.2	300–450	Subsoiler (straight shank-sweep tine, curved shank-chisel tine, curved shank-sweep tine, bentleg-chisel tine)	n/a	1230–1420	n/a	n/a	not stated
Kim et al. [44]	5.16 to 9.9	n/a	7.64–7.9	5–200	Mouldboard plough	Loam	1496–1904	24.5–34.02	n/a	EEPA
Aikins et al. [41]	5 to 31	8, 20 and greater	8	100	Bentleg and narrow point openers	Clay (Vertosol)	1504	23.7	46.4	HSCM + LCM
Wang et al. [74]	15.08	n/a	3	300	Winged subsoiler	Sandy loam	1404–1833	n/a	n/a	PBM
Sadek et al. [75]	≤20.2	n/a	4–16	127	Disc	Sandy loam	1700	16.32	n/a	PBM
Saunders et al. [76]	n/a	n/a	4.5–10	25–100	Plough skimmers	Sandy loam	1523.8	8.3	n/a	HSCM + LCM
Ma et al. [77]	2.88 to 5.97	n/a	1.08–2.16	120	Scraper	Sandy loam	1389	10	n/a	not stated
Hoseinian et al. [35]	2	2.5	0.9	150	Dual sideway-share	Sandy clay loam	1565	11.5	15.4	PBM

Table 2. Relative errors in DEM-predicted soil–tool reaction forces, travel speeds, and operating depths reported in the various literature across soil types.

Reference	Relative Error in DEM Prediction (%) Relative to Measured Data		Travel Speed (km h−1)	Operating Depth (mm)	Tillage Tools	Soil Texture	Dry Bulk Density (kg m−3)	Soil Water Content (%, w/w)	Cohesive Strength (kPa)	Contact Model
	Draught	Vertical Force
Sadek et al. [58]	n/a	n/a	n/a	n/a	n/a	Sandy soil	990 1280 1360 1500	0.02 13 21.5	1.23-32.70	PBM
Chen et al. [7]	4 to 31	n/a	3.19 (average)	100	Sweep tine	Coarse sand Loamy sand Sandy loam	1410 1330 1410	8.98 14.84 18.2	15.7 25.2 36	PBM
Obermayr et al. [52]	n/a	n/a	2.16–4.5	10–200	Bulldozer blade	n/a	1900	n/a	11.16	LSCM + cohesion
Tamas et al. [30]	4 to 12	n/a	1.8–8.64	200	Sweep tine	Sandy soil	1850	6.33	11.86	PBM
Bravo et al. [18]	9, 24	n/a	-	150–500	Para-plough and mouldboard plough	Clay (Vertosol)	1000 1200 1400	8 18 20 35	25–125	LSCM + cohesion
Li et al. [56]	3 to 15	n/a	3.6	180–260	Subsoiler	n/a	n/a	19	n/a	PBM
Mak and Chen [61]	n/a	n/a	2.2–6.59	50–200	Sweep tine	Loamy sand	1320	11.3	13.9	PBM
Obermayr et al. [72]	n/a	n/a		100–200	Straight-vertical blade and bulldozer blade	Sand	1520 1980 1870	10 15	6–22.5	LSCM + cohesion
Ucgul et al. [38]	≤11.6	≤15.2	5–12.5	70	Sweep tine	Sandy loam	1750	8	6	HSCM
Ucgul et al. [53]	n/a	n/a	4–12	75	Sweep tine	Sandy loam	1320 1780 1880	1 15 13	3 15 22	HSCM + LCM
Kotrocz et al. [60]	n/a	n/a	n/a	50–150	Cone penetrometer	Loamy sand	1632	15.8	6.61–8.66	PBM
Li et al. [70]			2.99	3–18	Claw	Sandy loam	1300	n/a	17.5	PBM
Murray [69]	1.86	50.7	8	38	Disc and hoe openers	Clayey lacustrine	1560	19.6	n/a	PBM
Hang et al. [32]	n/a	n/a	3	300	Subsoiler	Loamy clay	1346	12.5	11.8	PBM
Milkevych et al. [62]	n/a	n/a	3.2	100	Sweep tine	Coarse sand Loamy sand	1410 1330	9 14.8	15.8 25.1	PBM
Tekeste et al. [55]	9, 12	-59, -49	0.79–9.65	102	Sweep tine	Loam	1307	8.99	33	PBM
Tong et al. [73]	<10	<10	7.2	300–450	Subsoiler (straight shank-sweep tine, curved shank-chisel tine, curved shank-sweep tine, bentleg-chisel tine)	n/a	1230–1420	n/a	n/a	not stated
Kim et al. [44]	5.16 to 9.9	n/a	7.64–7.9	5–200	Mouldboard plough	Loam	1496–1904	24.5–34.02	n/a	EEPA
Aikins et al. [41]	5 to 31	8, 20 and greater	8	100	Bentleg and narrow point openers	Clay (Vertosol)	1504	23.7	46.4	HSCM + LCM
Wang et al. [74]	15.08	n/a	3	300	Winged subsoiler	Sandy loam	1404–1833	n/a	n/a	PBM
Sadek et al. [75]	≤20.2	n/a	4–16	127	Disc	Sandy loam	1700	16.32	n/a	PBM
Saunders et al. [76]	n/a	n/a	4.5–10	25–100	Plough skimmers	Sandy loam	1523.8	8.3	n/a	HSCM + LCM
Ma et al. [77]	2.88 to 5.97	n/a	1.08–2.16	120	Scraper	Sandy loam	1389	10	n/a	not stated
Hoseinian et al. [35]	2	2.5	0.9	150	Dual sideway-share	Sandy clay loam	1565	11.5	15.4	PBM

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 non-linear Hertzian spring model (Figure 6b) [79]. Figure 6 shows “a schematic diagram of particle contact and normal force-overlap (f_n-δ) 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,45,46].

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 time-consuming 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].

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,88,89,90,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,93,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.

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 m_s, the mass of the ball is m_r, and the angles at which sliding and rolling occur are Ψ_s and Ψ_r, respectively, then the coefficients are calculated according to Equations (28) and (29) [29,97].

Soil–tool coefficient of static friction,

$${\mu }_s=\frac{m_{s}g\sin {\Psi }_s}{m_{s}g\cos {\Psi }_s}=\tan {\Psi }_s$$

(28)

Soil–tool coefficient of rolling friction,

$${\mu }_r=\frac{m_{r}g\sin {\Psi }_r}{m_{r}g\cos {\Psi }_r}=\tan {\Psi }_r$$

(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 (F_a) is applied while an increasing horizontal (shearing) force (F_b) 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.

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 ϕ

(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 (σ₁) and radial (σ₃) stresses as indicated in Figure 12 until soil failure is achieved [98]. The radial stress (σ₃) 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 σ₃ to form a resultant axial stress σ₁. The above steps are repeated several times under increasing radial stress. The plots of the deviator stress (q = σ₁ – σ₃) 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].

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.

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.7^TM), 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 V_g being grid volume and V_p the total volume of particles whose centroids are located within the grid, voidage can be calculated according to Equation (31).

$$Voidage=\frac{V_g-V_p}{V_g}\times 100\%$$

(31)

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:

• 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).
• 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.

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].

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,107,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 (Figure 15 and Figure 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⁻¹ 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⁻¹ and 300 mm, respectively, had an RE of 9.71%, whereas the non-winged 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⁻¹ 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⁻¹).

6. Soils Modelled in DEM Simulations

Table 1 and Table 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⁻³) 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:

• Even though the Hertz–Mindlin contact model (HMCM) has been used in most DEM studies of tillage and furrow opening, it consistently fails to predict vertical soil force accurately. The Hysteretic Spring contact model (HSCM) can more accurately predict soil forces and particle movement.
• Angle of repose, inclined plane, direct shear, triaxial compression, and some in situ tests (grouser shear, plate sinkage, and cone penetration tests) have been used to measure and calibrate DEM input parameters. The angle of repose test has been used mainly for cohesionless soils due to the poor flowability of cohesive soils. However, using results from reproducible phases of the angle of repose experiment, successful calibrations for cohesive soils have been achieved.
• Unlike other numerical models, DEM is able to closely predict not only soil forces, but it is also capable of modelling soil failure mechanisms, soil loosening, and soil particle movement. Soil rupture and crack propagation, critical depth, three-dimensional particle movement within the soil profile and lateral particle movement on top of the soil have all been predicted in DEM.
• Using voidage or porosity grids to determine loosened furrow cross-sectional profiles has been found to be superior to using particle velocity and displacement profiles. However, some researchers have successfully used a particle displacement approach to determine accurate furrow profiles with a more objective criteria for defining loosened furrow boundary.
• Close predictions of draught and vertical forces (≤20%) have been obtained with DEM. These predictions can be improved by using smaller particles of a near-real shape. However, this must be balanced with computation time requirements.

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

• The Edinburgh elasto-plastic adhesion model (EEPA) has been successfully used to model consolidated or cohesive powders. This contact model is recommended to be studied more extensively for cohesive soils, although some researchers have used it.
• Due to pore water pressure within wet and soft soils, coupling DEM and CFD is likely to produce more accurate simulations. This idea can be explored in future research.
• A comprehensive analysis of soil disturbance parameters has been successfully done using voidage grids in EDEM^® DEM software. Replication of this approach in other DEM software is recommended.
• The criteria introduced by Aikins et al. [41] for defining particle displacement threshold for DEM furrow profile identification need further investigation with particles of smaller radii than the 5 mm used in the study. This approach can provide greater details on the three-dimensional soil translocation process.