Applications of Discrete Element Method in the Research of Agricultural Machinery: A Review

Abstract

As a promising and convenient numerical calculation approach, the discrete element method (DEM) has been increasingly adopted in the research of agricultural machinery. DEM is capable of monitoring and recording the dynamic and mechanical behavior of agricultural materials in the operational process of agricultural machinery, from both a macro-perspective and micro-perspective; which has been a tremendous help for the design and optimization of agricultural machines and their components. This paper reviewed the application research status of DEM in two aspects: First is the DEM model establishment of common agricultural materials such as soil, crop seed, and straw, etc. The other is the simulation of typical operational processes of agricultural machines or their components, such as rotary tillage, subsoiling, soil compaction, furrow opening, seed and fertilizer metering, crop harvesting, and so on. Finally, we evaluate the development prospects of the application of research on the DEM in agricultural machinery, and look forward to promoting its application in the field of the optimization and design of agricultural machinery.

1. Introduction

In order to improve the operation performance and efficiency of agricultural machines and satisfy the production requirements of modern agriculture, it is necessary to upgrade, update, design, and optimize agricultural machines and their components [1,2]. Theoretical analysis, numerical calculation, and experiment are the three most commonly used approaches in the design and optimization process [3]. Computer simulation, which shows the great advantages of a short cycle, low cost, and being free from the farming season compared with experiments, has been popularly used in recent years in the research of agricultural machinery [4,5].

There are two generally adopted computer simulation methods, the finite element method (FEM) and the discrete element method (DEM). FEM has been mostly used in the research of continuous mediums but most agricultural materials are discrete, such as soil particles and crop seeds; therefore, FEM shows limitations in the study of soil flow and mixture, failure and deformation, and the motion of seed flow [6,7]. Whereas DEM, a numerical calculation approach for the analysis of complex dynamic discontinuous mechanical discrete systems [8,9], has shown its advantages in the research of agricultural machinery.

In DEM, research targets are divided into the discrete elements of particles and walls. In the simulation of operational process of agricultural machinery, agricultural materials can be modeled by particles or particle agglomerates, while agricultural machinery or its components can be modeled by a combination of walls. The trajectory, movement, and mechanical behavior of the agricultural materials and the resistance of the agricultural machinery and its components can be recorded in real time, so as to conduct an in-depth study of the interaction mechanism between agricultural materials and agricultural machinery, and offer guidance for the design and optimization of agricultural machinery.

Recently, many kinds of DEM software have been developed, such as EDEM, PFC, Pasimodo, Agridem, Yade-DEM, DEMeter++, and so on [10,11,12,13]. Using this software, many kinds of agricultural material models with high accuracy have been established, and with these established models, many simulation studies of the interaction processes between agricultural machinery and agricultural materials with high reliability and high fidelity have been carried out [14,15].

The aim of this review is to introduce the application research progress of DEM in agricultural machinery, including the model establishment of common agricultural materials and the typical operational process simulation of agricultural machinery and its components. The prospects of DEM’s application in research on agricultural machinery are then proposed, while looking forward to promoting the application of DEM in the field of the optimization and design of agricultural machinery.

2. DEM Model Establishment and Parameter Calibration of Agricultural Materials

Accurate and reliable DEM models of agricultural materials are the foundation of the high fidelity simulation of agricultural machinery operation processes. Agricultural materials mainly include soil, crop seed, straw, and fertilizer, etc.

2.1. Soil Model Establishment and Parameter Calibration

Soil is one of the basic research objects in agricultural machinery. Affected by geographical location, climate, soil texture, and many other factors, the soil in different regions shows a great variance in physical and mechanical properties. As listed in Table 1, various contact models have been developed in DEM to represent the contacts between soil particles with different physical and mechanical properties, providing the basis for the establishment of a soil model. To establish a soil model, a suitable contact model needs to be established first, then the key parameters that affect the mechanical and dynamic behavior of the soil need to be calibrated.

2.1.1. Establishment of Soil Contact Model

A DEM soil model can be classified into a cohesive model or a cohesionless model, based on whether the soil has cohesion [16]. In the modeling of cohesionless soil, the Hertz–Mindlin contact model (HMCM), in which the deformation at the contact point is assumed to be non-linear elastic, can simulate the movement of cohesionless soil and the resistance of the agricultural machinery with a certain accuracy, but only the elasticity behavior of the soil can be modeled [17,18,19]. Ucgul et al. [20] developed a hysteretic spring contact model (HSCM), consisting of linear elasticity and cohesion, where the contact between particles of this model showed a linear elasticity when the stress was low, once the stress exceeded the threshold value, particle contact behaved in a plastic way. Experiment results showed that compared with HMCM, HSCM improved modeling accuracy in operation resistance under different travelling speeds, especially in vertical resistance [21]. Moreover, in order to simplify the calculation process, the contact between soil and machine mainly adopted a cohesionless contact model such as the Hertz–Mindlin (no slip) model [22,23].

Soil shows cohesion due to the existing liquid bridge. The parallel bond contact model (PBCM), which takes cohesion into consideration based on HMCM, is often adopted in DEM software PFC (Itasca Consulting Group, Inc., Minneapolis, MN, USA) [24,25]. The cohesion between soil particles was represented by a cylinder material, through which both force and moment can be delivered. When the distance between two adjacent particles was smaller than the preset value, the parallel bond was created spontaneously; when the stress exceeded a threshold value, the parallel bond broke off, and the two particles were separated, so as to mimic the soil breaking behavior. PBCM has mostly been adopted in sandy soil or sandy loam soil models, but only cohesions between soil particles can be simulated [26,27]. By adding cohesion and adhesion in a normal direction in the HSCM, the adhesion between soil and machinery can be modeled [28]. Experiments showed that the improved HSCM with LCM (liner cohesion model) could predict resistance more accurately under different moisture contents and pressures [29,30]. The Hertz–Mindlin with bonding model, one of the most used models in EDEM software (DEM Solutions Ltd., Edinburgh, UK), employed an “adhesive glue” to stick soil particles together. Shear and normal displacement could be carried with the “adhesive glue” until the shear stress reached the maximum value, then the “adhesive glue” breaks, in this way soil breaking behavior in rotary tillage, subsoiling, and separation of soil and potatoes and so on could be simulated [31,32,33]. However, this model was only appropriate for soils with low cohesion, because the soils with a high cohesion show more cohesive and plastic behavior in a macroscopic view [34]. The Hertz–Mindlin with JKR cohesion model uses JKR (Johnson–Kendall–Roberts) normal elasticity contact force to represent soil cohesions [35]; with this model, Xing et al. [36] established a model of lateritic soil in hot areas of Hainan Province. The soil breaking resistance error was 3.43%; Li et al. [34] simulated a clayey black soil in the Northeast Plain in China, and the soil aggregation phenomenon when moisture content exceeded 12% could be simulated.

All the above-mentioned soil models were constituted by spherical particles with the same physical and mechanical properties. However, soil composition is complex in real fields. Soil properties show distinctions at different soil depth, and are often mixed with straw, stubble, and other materials. Current studies have established a compound model of different soil layers and a compound root–soil model. In the establishment of a compound soil model of different soil layers, one method is to set different colors as different soil depth layers (Figure 1a) [37,38,39]; the other is determining the soil contact and constitutive parameters according to the soil structure, such as the tillage layer, plow-pan layer, and subsoil layer for subsoiling (Figure 1b) [23]. These models can be used to simulate the movement of soil particles at different depths and the resistance produced by specific soil layers on the tillage device. In the study of a root–soil model, Liu et al. [40] established a maize root model with bonded DEM particles and fixed the root model with soil particles. The model could be broken off by external forces (Figure 1c); the root–soil model established by Frank et al. [41] used sphere particles to model soil and flexible cylinders to model stubble, the stubble could bend and stretch, and the soil could be broken off. This root–soil model could be used in the study of tillage processes in a no-till field. These compound models offer new concepts for the modelling of complex soil.

2.1.2. Calibration of Physical, Contact, and Constitutive Parameters of a Soil Model

In order to predict the movement of soil particles and resistance in the operation of agricultural machinery accurately, the DEM parameters of a soil model should be calibrated.

There are two kinds of DEM parameter, namely physical parameters, such as porosity and density, and contact and constitutive parameters, such as the friction coefficient, particle strength and stiffness, and bond strength and stiffness [19]. Most research has modeled soil particles with sphere particles, while Tekeste et al. [42] scanned the soil profile with a 3D scanner (Artec Space Spider model, Artec Studio, Luxembourg), then imported the profile data into a 3D modeling software ANSYS SpaceClaim (SpaceClaim Corp., Concord, MA, USA) and combined spherical particles to create a soil particle model that resembled the real profile. The radius of soil particles has a vital influence on the modeling accuracy and time spent. Many studies found that a decreased particle radius could improve modeling accuracy, with an increase in the time spent [21,28,43]. For the purpose of reducing the time spent, it was common to increase particle radius, so long as it could still ensure modeling accuracy [44,45,46,47]. When the particle radius is changed, other parameters need adjustment at the same time. Ucgul et al. [20] found that when soil particle radius was increased from 4 to 9 mm, the collision coefficient of restitution, static friction coefficient, rolling friction coefficient, and time step all increased; Thakur et al. [48] found that in order to match the macroscopic motion behavior with real experiments, particle stiffness should be increased after the increase of particle radius.

There are two common methods for calibrating the contact and constitutive parameters [15,49]. The first is determining the parameters directly with test apparatus; the second is comparing the simulated dynamic or mechanistic behaviors with the experimental ones through a repose angle test (Figure 2a), direct shear test (Figure 2b), triaxial pressure test (Figure 2c), and tillage test, by changing the parameter value until the error between the simulation and experimental results was controlled in a limited range, optimal parameter values could be obtained [17,27,29,50,51,52].

Friction coefficient is one of the most important parameters influencing soil dynamic behavior, including the static and rolling friction coefficient between soil particles, and the static and rolling friction coefficient between soil particles and the contact materials. The repose angle test was the most commonly used method to calibrate friction coefficient between soil particles. First, a repose angle experiment was conducted to obtain the real repose angle; second, the friction coefficient value of the DEM model was changed until the simulated repose angle was approximately that of the experiment value, then the friction coefficients could be calibrated [53,54]. Friction coefficients between soil particles and the contact material could be obtained directly with a sliding test on a slope or flat sliding test [53]. The collision coefficients of the restitution between soil particles and contact material were generally determined by collision test [55]. Soil cohesion was mostly affected by particle stiffness, particle strength, bond stiffness, and bond strength [56,57,58,59,60]. Mak et al. [22] found that the operation resistance of a simple soil engaging tool increased with the increase of particle normal and shear stiffness, the soil disturbance area increased with bond shear strength, and the optimal DEM parameter values varied with soil type. Chen et al. [53] calibrated particle stiffness for three kinds of soil by comparing the draft force, vertical resistance, and soil disturbance in a furrow opening test; the relative error was smaller than 10%. Pue et al. [58] calibrated the soil Young’s modulus, Poisson’s ratio, and bond shear and normal stiffness with a triaxial compression test.

Table 1 summarizes the commonly used soil contact models and key model parameters in the published literature.

Table 1. Common DEM models and key parameters of soil.

Soil Contact Model	Soil Type	Key DEM * Parameters
		Static Friction Coefficient between Soil Particles	Rolling Friction Coefficient between Soil Particles	Static friction Coefficient between Soil Particles and Contact Material	Rolling Friction Coefficient between Soil Particles and Contact Material	Collision Coefficient of Restitution	Mechanical Parameters
HMCM *	sea sandy soil [16]	[18]			Steel [53]		Elasticity modulus [18] Particle stiffness [16]
Hertz-Mindlin with bonding	Norfolk sandy loam soil [52] Clay soil in South China [31] Loam clay soil [32]	[53,54]	[53]	Steel [53] Steel [54]	Steel [54]	Between soil particles [55] Between soil particle and contact material [55]	Shear modulus [53] Poisson’s ratio [53] Bond normal stiffness [52] Bond shear stiffness [52] Critical stress of the bond [32]
PBCM *	Sandy loam soil [22,56] Sandy soil [56] Sandy soil [27,57]	[56]	[56]				Normal and shear particle stiffness [22,56] Cohesion strength [58] Young modulus [27] Elasticity modulus of particle; elasticity Modulus of bond; damping coefficient [59]
HSCM *	Sandy loam soil [28]	[28]	[28]	Steel [28]	Steel [28]
HSCM with LCM *	Sandy loam soil in North-west area [29] Grapevine antifreezing soil in Ningxia, China [30] Light loam in North China Plain [23]	[23,29]	[23,29]	Steel [23]	Q235 steel rubber [30] Steel [23]	Between soil particle and Q235steel/rubber [30]	Shear strength, yield strength [29] Surface energy [23]
Hertz-Mindlin with JKR * Cohesion	Clayey black soil [34] Lateritic soil particles in hot areas [36]	[34]	[34,36]	65 Mn; PTFE * [34]	65 Mn; PTFE [34] 28 Mn; B5 board [36]	Between soil particles and 65 Mn/PTFE [34] Between soil particles [36]	Surface energy between soil particles and 65 Mn/PTFE [34] Surface energy between soil particles [36] Poisson’s ratio [34]

* Note: Define of the abbreviations. DEM (discrete element method), HMCM (Hertz–Mindlin contact model), PBCM (parallel bond contact model), HSCM (hysteretic spring contact model), LCM (liner cohesion model), JKR (Johnson–Kendall–Roberts), PTEE (polytetrafluoroethylene).

Table 1. Common DEM models and key parameters of soil.

Soil Contact Model	Soil Type	Key DEM * Parameters
		Static Friction Coefficient between Soil Particles	Rolling Friction Coefficient between Soil Particles	Static friction Coefficient between Soil Particles and Contact Material	Rolling Friction Coefficient between Soil Particles and Contact Material	Collision Coefficient of Restitution	Mechanical Parameters
HMCM *	sea sandy soil [16]	[18]			Steel [53]		Elasticity modulus [18] Particle stiffness [16]
Hertz-Mindlin with bonding	Norfolk sandy loam soil [52] Clay soil in South China [31] Loam clay soil [32]	[53,54]	[53]	Steel [53] Steel [54]	Steel [54]	Between soil particles [55] Between soil particle and contact material [55]	Shear modulus [53] Poisson’s ratio [53] Bond normal stiffness [52] Bond shear stiffness [52] Critical stress of the bond [32]
PBCM *	Sandy loam soil [22,56] Sandy soil [56] Sandy soil [27,57]	[56]	[56]				Normal and shear particle stiffness [22,56] Cohesion strength [58] Young modulus [27] Elasticity modulus of particle; elasticity Modulus of bond; damping coefficient [59]
HSCM *	Sandy loam soil [28]	[28]	[28]	Steel [28]	Steel [28]
HSCM with LCM *	Sandy loam soil in North-west area [29] Grapevine antifreezing soil in Ningxia, China [30] Light loam in North China Plain [23]	[23,29]	[23,29]	Steel [23]	Q235 steel rubber [30] Steel [23]	Between soil particle and Q235steel/rubber [30]	Shear strength, yield strength [29] Surface energy [23]
Hertz-Mindlin with JKR * Cohesion	Clayey black soil [34] Lateritic soil particles in hot areas [36]	[34]	[34,36]	65 Mn; PTFE * [34]	65 Mn; PTFE [34] 28 Mn; B5 board [36]	Between soil particles and 65 Mn/PTFE [34] Between soil particles [36]	Surface energy between soil particles and 65 Mn/PTFE [34] Surface energy between soil particles [36] Poisson’s ratio [34]

* Note: Define of the abbreviations. DEM (discrete element method), HMCM (Hertz–Mindlin contact model), PBCM (parallel bond contact model), HSCM (hysteretic spring contact model), LCM (liner cohesion model), JKR (Johnson–Kendall–Roberts), PTEE (polytetrafluoroethylene).

2.2. DEM Model Establishment of Crop Seed and Parameter Calibration

The profile of crop seeds is usually irregular, the main process of model establishment is: determine the profile and establish a 3D outline model; establish a DEM seed model with particles; calibrate the key parameters affecting seed dynamic and mechanistic behaviors and verify the results with experiments. Common crop seeds have no cohesion between each other in the operation of agricultural machinery; a Hertz–Mindlin (no slip) model, which has no cohesion, has been widely adopted [4].

2.2.1. 3D Model Establishment of Crop Seeds

There were three common approaches to determining crop seed profile: direct determining, slice modeling, and 3D-scanning. The direct determining method means determine the length, width, and thickness of crop seeds with a micrometer or vernier caliper directly, and calculating the equivalent diameter and degree of sphericity [61,62,63]. Based on the difference of seed profile and size, a model can be established by classification [64,65,66]. The slice modeling method slices crop seeds into a certain thickness and collects an outline profile of each section, then imports them into 3D modeling software to establish a 3D model [67,68]. The 3D-scanning method uses 3D laser scanner or minitype CT scanner to scan the seed profile and obtain the point cloud data, then imports the data into 3D modeling software. This was more accurate than the other two approaches [69,70]. By these methods, crop seed models of wheat, maize, rice, potato, flax, oilseed rape, soybean, and black pepper have been established in recent studies, as listed in Table 2 [64,65,67,71,72,73,74,75,76,77,78,79].

A crop seed model is usually established by a particle agglomeration method, and most studies found that the simulation accuracy increases with the increase of particle quantity of the model, but the time spent increased simultaneously [69,80]; however, some researchers found that specific crop seed models were more accurate when the particle number was less [71,81].

2.2.2. The Contact and Constitutive Parameters of a Crop Seed Model

Apart from outline profile, the DEM parameters of a crop seed model also include contact parameters such as friction coefficient, elasticity modulus, and so on. The repose angle test was widely used to determine the static and rolling friction coefficient between crop seeds, and its calibration process was similar with the soil parameters. As listed in Table 1, the static and rolling friction coefficients of many kinds of crop seed have been calibrated [78,82,83,84,85,86,87,88,89,90,91,92]. The static and rolling friction coefficients between crop seeds and contact materials were usually determined with a slope sliding method or self-made coefficient determining apparatus [64]; collision coefficients of restitution were calibrated with seed dropping tests and collision tests [74,75]. Using a static compression test, stiffness coefficient, elasticity modulus, Poisson’s ratio, and shear modulus were calibrated [62,72,93]. Some other studies calibrated the damping coefficient with a uniaxial loading test [94,95,96]. Common crop seed models and the key DEM parameters are listed in Table 2.

In real operation processes, air resistance and magnetic forces may also exist. Binelo et al. [77] established an aerodynamic engine, through which the seed suspending velocity, drag coefficient, orientation, and projection area were added in the model. In this way, the air friction of the total seed could be calculated. In order to model the magnetic force of seeds coated with magnetic powder, Wang [97] and Hu et al. [98] created a force plug-in with C programming language using numerical fitting, then loaded the magnetic power through an API interface of EDEM software (DEM Solutions Ltd., Edinburgh, UK) to model the magnetic attraction behavior. Yu et al. [99,100,101] established a maize-ear model with Agridem software, which they developed. A corncob was modeled by particles with connections, and special balls were set between the niblet and corncob particle to connect them together, the connections could be broken off, so that the breakage of the corncob and maize threshing process could be simulated.

Table 2. DEM model and key parameters of common crop seeds.

Seed Type	3D Model	Key DEM Parameters
		Coefficient between Soil Particles			Coefficient between Soil Particles and Contact Material			Elasticity Modulus	Damping Coefficient
		Static Friction	Rolling Friction	Collision Restitution	Static Friction	Rolling Friction	Collision Restitution
Wheat	[71,72]	[71,82]	[71]	[82]	Steel [82] Low-carbon steel [71]	Low-carbon steel [71]	[82]		Between wheat seed; wheat seed and steel [82]
Maize	[64,65,66,67]	[64,83]	[64,71,83]		Crylic acid [83] Low-carbon steel [71]	Crylic acid [83] Aluminium; Organic glass [90] Low-carbon steel [71]
Rice	[73] [74]	[74]	[74]				[74]
Potato seed	Sphere; big ellipsoid particle; small ellipsoid [65]	[84]		[65]	Steel; plastic [65]	Steel; plastic [84]	Steel [65]
Flax	[75]	[75]	[75]	[75]	Organic glass; Aluminum cylinder [75]	Organic glass; Aluminum cylinder [75]	Organic glass; steel [75]	[75]
Rapeseed	[76]	[85,86]	[85,86]	[85]	Abs; aluminium alloy [85] Organic glass [87]	Abs; aluminium alloy [85] Organic glass [76,87]	Abs; aluminium alloy [85]	[95]	[95]
Soybean	[77,78]	[78]	[78]				[78]
Black pepper	[79]	[79]	[79]		steel [79]	steel [79]	[79]

Table 2. DEM model and key parameters of common crop seeds.

Seed Type	3D Model	Key DEM Parameters
		Coefficient between Soil Particles			Coefficient between Soil Particles and Contact Material			Elasticity Modulus	Damping Coefficient
		Static Friction	Rolling Friction	Collision Restitution	Static Friction	Rolling Friction	Collision Restitution
Wheat	[71,72]	[71,82]	[71]	[82]	Steel [82] Low-carbon steel [71]	Low-carbon steel [71]	[82]		Between wheat seed; wheat seed and steel [82]
Maize	[64,65,66,67]	[64,83]	[64,71,83]		Crylic acid [83] Low-carbon steel [71]	Crylic acid [83] Aluminium; Organic glass [90] Low-carbon steel [71]
Rice	[73] [74]	[74]	[74]				[74]
Potato seed	Sphere; big ellipsoid particle; small ellipsoid [65]	[84]		[65]	Steel; plastic [65]	Steel; plastic [84]	Steel [65]
Flax	[75]	[75]	[75]	[75]	Organic glass; Aluminum cylinder [75]	Organic glass; Aluminum cylinder [75]	Organic glass; steel [75]	[75]
Rapeseed	[76]	[85,86]	[85,86]	[85]	Abs; aluminium alloy [85] Organic glass [87]	Abs; aluminium alloy [85] Organic glass [76,87]	Abs; aluminium alloy [85]	[95]	[95]
Soybean	[77,78]	[78]	[78]				[78]
Black pepper	[79]	[79]	[79]		steel [79]	steel [79]	[79]

2.3. Other Agricultural Materials

The establishment of a straw model has been the latest development in DEM research. Straws physical structure is complicated; in the preliminary stage, a rigid straw model was mostly studied, which could be used to analyze the straw movement, burying, and distribution in the operation of soil tillage [102,103,104,105,106]. By adding contact points or a stick key between two rigid straw models, moment and force could be transferred, and a flexible straw model with elasticity was established [107,108,109]. Furthermore, Schramm et al. [110], Wang et al. [111], and Liu et al. [112] calibrated the Young’s modulus, cohesion damping coefficient, elasticity modulus, and cohesion parameters through a cantilever beam test and three point bending test. For the straw cutting behavior, Guo [113] and Zhang [114] modeled the cutting process of rattan straw and maize straw scarfskin, respectively. Liao et al. [115,116] established fodder rape straw models in the bolting stage and early pod stage, respectively, and calibrated the friction coefficient and cohesion parameters through a repose angle test and straw cutting test. Using a straw compression test and cutting test, Zhang et al. [117] calibrated the key mechanic parameters of Young’s modulus, bending strength, and elasticity modulus of a BPM (bonded particle model) straw model; the model could simulate the four straw forms of short, standard, long, and unbroken straw after being crushed [118,119,120].

Table 3. DEM models and key parameters of straw.

Straw Type	Straw Model	Specialty	Contact Parameters	Constitutive Parameters
Oat straw	[104]	Rigid	Friction coefficient between straw and blade [104]	Straw stiffness [104]
Maize straw	[117]	Cuttable	Normal and shear contact stiffness, in normal and shear critical stress [117]
Fodder rape straw in bolting stage	Simulation Experiment [115]		Collision coefficient of restitution between straw and between straw and steel; static and rolling coefficient between straw and between straw and steel [115]	Normal and shear stiffness; normal and shear critical stress; bond radius [115]
Fodder rape straw in early pod stage	[116]		Collision coefficient of restitution between straw and between straw and steel; static and rolling coefficient between straw and between straw and steel [116]	Elasticity modulus, shearing modulus and Poisson’s ratio [116]
Wheat straw	[110]	Flexible	Viscous damping coefficient [110]	Young’s modulus [110]
	Experiment Simulation [109]		Friction coefficient between straw, Friction coefficient between straw and steel [118]	Young’s modulus [118,120] Bending strength [108] Tension modulus [120]
	[112]		Static friction coefficient between straw and steel, restitution coefficient [112]	Shearing modulus; elasticity modulus Bond radius, shear and normal cohesion stiffness [112]

summarizes the established straw models and their key DEM parameters in the literature.

Fertilizer models were mostly represented by single spherical particles: the static and rolling friction coefficient between urea particles and the static and rolling friction coefficient between urea particles and ABS (acrylonitrile-butadiene-styrene) and PVC (polyvinyl chloride) material were calibrated through repose angle test and slope sliding test [121,122]; and the terminal velocity of large particle urea, DAP (diammonium phosphate), and potassium sulfate was determined through CFD-DEM coupling simulation [123]. Based on the mix uniformity of three fertilizers after fertilizing, Yuan et al. [124] calibrated the rolling friction coefficient between nitrogen fertilizer, phosphate fertilizer, and potassic fertilizer, and calibrated the static and rolling friction coefficient between the fertilizer and a conveyor.

There have also been studies that established a DEM model of fodder, vermicomposting nursery substrate, and pig manure organic fertilizer treated with Hermetia illucen, and calibrated the contact and constitutive DEM parameters [125,126,127]. Coetzee and Lombard [128] established a grape model to predict the removal of berries from the stems (Figure 3); grape berries were modeled by single DEM particles, and the stem was modeled by a branch of particles that stick together and which could be broken off, and the model could simulate the grape picking and stem breaking process.

3. DEM Simulation of Agricultural Machinery Operation Processes

The interaction relationship between agricultural material and agricultural machinery during its operational process has been the focus of DEM modeling research [14]. The DEM simulations of typical agricultural machinery operation processes are as follows:

3.1. Simulation of Tillage and Soil Preparation

3.1.1. Simulation of Rotary Tillage

DEM could record the resistance and torque of a rotary blade in a small time interval, the variation trends of the simulated working resistance curve in all directions were drawn, and they were found to be similar to the theoretical resistance curve and the soil bin experimental results [129,130,131]. For torque requirements, the change of key stages can be obtained and typical variation curves can be illustrated with DEM (Figure 4) [130,132,133,134]. The experiment curve had a second peak value, while the simulation one did not, indicating the simulation remains to be improved [135].

For the aspect of soil and straw disturbance, DEM software can monitor the path, contact force, displacement, velocity, and accelerated velocity of single particles, and analyze their dynamic characteristics under different working depths and rotating speeds microscopically [31,102] (Figure 5a). It can also record the soil and straw disturbance produced by a rotary tillage device quantitatively in a macroscopic view [60,103], as shown in Figure 5c,d; Zheng et al. [136] marked soil blocks in different colors in a DEM model and found that a rotary blade could shear, tear, overturn, and throw the soil block so as to make the tillage layer flat, comminute, and loosen. Zhao [137] compared velocity vector maps of different types of rotary blade (Figure 5d) and found that when the blade was tilted with the travel direction, straw could be thrown aside so that the seedbed could be cleaned up.

3.1.2. Simulation of Subsoiling

The purpose of subsoiling is to break up the plow-pan and loosen the soil. Zeng et al. [3] used a built-in measuring sphere in DEM to monitor the internal stress and porosity, with which the break up performance of a plow-pan was evaluated; while Ding et al. [32] evaluated the soil breakage performance by calculating the number of broken bonding keys of a DEM cohesive soil model (Figure 6a). Huang et al. [138] analyzed the microcosmic movement and macroscopic disturbance of soil particles in different positions, and clarified the soil dynamic principles during a subsoiling process (Figure 6b). By combining a DEM simulation with an orthogonal test method, the effect of the structural parameters of the subsoiler, such as if it has a wing, the install height of the wing, the distance between two subsoilers, and the operational parameters (subsoiling depth, travel speed on working resistance, and soil disturbance) have been studied [139,140,141,142]. Based on the above-mentioned results, various subsoilers, such as a polyline soil-breaking blade subsoiler and convex blade subsoiler (Figure 6c), have been developed [23,143,144,145,146].

3.1.3. Simulation of Plowing

The aim of plowing is to bury the surface soil and straw into the deeper soil layer. With the established compound soil model with different layers, Ucgul et al. [45,47] recorded soil particle movement at different depths. It was found that when the plowing depth was 300mm, only 10% of the surface soil remained at a 100 mm depth, 53% soil particles were buried at 200–300mm depth, and the simulation results were verified with field experiments. Saunders et al. [147] studied the relationship between the soil burying performance of a moldboard skimmer and the drag force, and found that DEM could predict working resistance more accurately than the other existing theoretical models.

3.1.4. Simulation of Interrow Tillage

In the simulation research of interrow tillage, by monitoring soil particle movement and the breakage of bonding keys, Cheng [148] analyzed the influence of rotary speed, travel speed, bending angle of the hiller blade, number of hiller blades, angle of the shell and soil breakage, number of soil particles being thrown upward, and distribution of the thrown soil. Liu [149] studied the influence of tillage depth, working speed, and angle of sweep wing on soil breakage rate, hilling height, and hilling width. This research was helpful for understanding the interrow tillage mechanism and provided reference for the design and optimization of structures and operation parameters of its key components.

3.1.5. Simulation of Soil Compaction and Traction

Soil is compacted during the operation of agricultural machinery. DEM could model the compaction of tires on soil, the sinkage of the tire, and the resistance and torque requirement of the tire [150,151]. Difficulties lie in the modeling of the pressure delivered from the tire to the soil surface, and three approaches are commonly adopted at present. The first was to model the pressure delivery with a built-in Sohne model (Figure 7a), with which the stress distribution could be simulated [152,153,154]. The second established a tire model with a particle agglomeration method, so that the tire has gravity and its sinkage could be simulated (Figure 7b) [155]. The third was with the help of a CFD (computational fluid dynamics)–DEM coupling method, where a tire was modeled with CFD and the soil was modeled with DEM (Figure 7c), and the dynamic characteristics of the tire on a sandy soil road could be simulated [156,157].

Traction performance can be simulated and evaluated by DEM. Nakanishi et al. [158] investigated the effect of lug height, lug thickness, number of lugs, and wheel diameter on the net traction of a wheel. Zeng et al. [159] predicted the tire rutting, vertical force, tractive force, and sinkage of a tire–sand interaction with DEM–FEM, and it was in good agreement with soil bin experiment. Nishiyama et al. [157] proposed a FEM–DEM coupling method for tire traction analysis, in which the soil model was transferred from a FEM to DEM model only when the tire was approaching, once the tire left, the model became a FEM model again. With this approach, the simulation time could be reduced significantly.

3.2. Simulation of Seeding

3.2.1. Simulation of Furrow Opening

With the help of DEM, recent studies have conducted single parameter tests to analyze the effect of working depth, rake angle, and travel speed on furrow opening performance and resistance [160,161,162]. By establishing opener models with different structures: Tamas et al. [53] compared the operation performance of different furrow openers; Tekeste et al. [40] studied the effect of worn openers with different surface profiles on working resistance; and Ucgul et al. [28,163] studied the effect of soil characteristics on furrow opening performance by changing the moisture content and compaction level of the soil model. Combined with an orthogonal test, Zhang et al. [164] and Liu et al. [165] optimized the key structure parameters of an opener and the moldboard shape of a furrow opening device, respectively.

3.2.2. Simulation of Seed Metering

In DEM simulation, the motion state of seed flow can be monitored in real time, and the effect of structure parameters and operation parameters of metering device on the metering performance can be evaluated [166] (Figure 8); for example, the effect of structure and operation parameters of a metering wheel on the beginning angle and ending angle of seed cleaning performance [167], and the motion states of maize seed groups being affected by vibration in field operations [168]. The force transfer of seed flow can also be monitored. Tian [169] studied the variation of the maximum force of rice seed in the seed taking stage of an ejection ear spoon device, which was difficult to obtain through an experiment method, and analyzed the reason for seed miss-filling and re-filling. Taking the metering percentage of a pass as the response target, a lot of structure and operation parameters have been designed and optimized with the help of DEM, such as depth, length, and section size of the model-holes of a hill-seeding centralized metering device [170]; and the arc radius, center angle, and lateral spacing of the middle plate of the seed taking claw of a garlic seed-picking device [171].

The CFD–DEM coupling method, which models airflow with CFD and models crop seeds with DEM, could monitor not only the airflow velocity and distribution in the metering device and metering tube, but also the movement of seeds in the airflow, so as to analyze: the effect of the structure of a pressurized tube on seed motion, and airflow in the tube of an air-assisted centralized metering device for rapeseed and wheat (Figure 9) [72,172]; the effect of the structure of three different metering discs of inside-filling air-blowing seed metering device on the drag force of the seed in a large circular shape [173]; the effect of the diameter and length of metering tube on airflow and motion of seed flow [174]; the seed filling performance in the adsorption, removing, and separation stages of a pneumatic seed metering service, with guided assistant filling [175], providing references for the design and optimization of air-assisted metering devices.

3.2.3. Simulation of Fertilizer Metering

In DEM, it is convenient to record the trajectory of fertilizer particles in the metering process, determine the spatial distribution position, distribution evenness index, etc. With the obtained data, a set of factor influencing the movement of fertilizer particles and fertilizer metering performance can be analyzed, such as the feeding rate, feeding angle and feeding position angle’s effect on the distribution in lateral distance of a fertilizer spreader with centrifugal swing disk [176] (Figure 10a); the effect of the inner diameter of the screw blade, screw pitch, outlet distance, number of screw heads, and the blocking wheel opening width on the coefficient of variation of fertilization stability [177].

With CFD–DEM coupling simulation, Gao [178] studied the effect of the rotary speed of a metering shaft on the metering quantity of a metering device, and analyzed the effect of position of the maximum airflow velocity and airflow velocity in the entrance on the movement of fertilizer in the tube; Liu et al. [179] analyzed airflow distribution and movement of a fertilizer group under different inlet velocities by monitoring the fertilizer accumulation performance through DEM, and evaluated the injection performance (Figure 10b).

Zeng et al. [5] monitored the dynamic characteristics between cotyledon and soil particles in soybean seedling emergence in DEM. Zhou et al. [180] analyzed the contact between seed and soil in a quantitative way, by monitoring their contact numbers, and providing new ideas for the understanding of seedbed preparation.

3.3. Simulation of Crop Harvesting

3.3.1. Simulation of Material Transfer

DEM modeling of the crop harvesting process includes material transfer, threshing, and cleaning. In order to analyze the wheat harvesting process, Wang et al. [181] adopted a coupling simulation method with EDEM–Recurdyn (EDEM, DEM Solutions Ltd., Edinburgh, UK; Recurdyn, FunctionBay Inc., Seoul, Korea), in which the wheat was simulated by DEM and the harvest combine was simulated by FEM, and the translate characters, velocity in an axial direction, and regional flow rate of material quantity in the continuous delivery process of the wheat were analyzed (Figure 11). Wang et al. [182] modeled rice plants with connected balls and simulated the transfer process from the conveyer to the threshing units of a combine harvester.

3.3.2. Simulation of Threshing

With the established corncob DEM model, Yu [100] modeled the maize threshing process of a drum-type corn threshing device (Figure 12), and analyzed the effect of rotating speed and feeding rate on the threshed rate and force of the corncob. Mou et al. [183] modeled the breaking process of maize kernels in a silage harvesting process, and obtained the optimal parameter combination of teeth number, blade depth and crushing gap, and the rotary speed of the knife roll of the threshing device. With the threshing simulation of rice, the large deformation and fragmentation of the rice plants was modeled by Wang et al. [182].

3.3.3. Simulation of Cleaning

To study the separation of straw and grain in the screening process, Lenaerts et al. [107] analyzed the properties of straw and grains on screening speed (Figure 13a). Li et al. [184] and Wang et al. [185] effectively analyzed the influence of the operation parameters of vibrational amplitude, frequency, and direction angle on screening time and efficiency. Han [186] and Ma et al. [187] optimized the structure of separation sieves (Figure 13b).

Using a CFD–DEM coupling method, the movement of particles in the air-screen cleaning was monitored, and the effect of inlet airflow velocity on the longitudinal velocity and vertical height of grains and short straws was investigated [188]. Xu et al. [189] studied the centroid velocities of grains, stems, and light impurities in the air-and-screen cleaning process, and monitored their degree of dispersion. With these data, the cleaning performance was analyzed.

Wei et al. [33] modeled the separation of soil and potatoes during the potato harvesting process and clarified the effect of structure and operation parameters of a wavy separating sieve on the collision between potatoes and soil clod breakage.

3.4. Simulation of Post-Harvest Processes

3.4.1. Simulation of Grain Conveying

With a DEM simulation, it was observed that the rotation angles of paddy grains increased with the increase of feeding direction. Adjusting the grain vertical direction into the hulling region by controlling the feeding angle could increase the one-time hulling ratio and reduce the breakage ratio (Figure 14) [190]. Chen et al. [191] also found that paddy rice, not only has a translational motion, but also rotates due to the shear of grain flow, which leads to a change in the orientation angle, which was determined by the depth of the grain layer. This study could provide guidance for the design and optimization of feeding systems.

Wang et al. [192] studied the effect of the inclination angle of returning plate and belt travel speed on particle return rate of wheat granules in a returning device for a fully enclosed belt conveyor; the returning mechanism was clarified and the optimal combination parameters of the two factors was obtained.

3.4.2. Simulation of Storage and Discharge of Seed in Silos

With DEM, the discharge process of a hopper was studied, it was found that the shape of the particle had a significant influence on discharge rate [193]. Using the sum of vertical forces exerted by rapeseed particles on the walls and floor of a receiving container, a difference in the mass flow rates was observed [87]. By recording the position, velocity, and force of wheat particles in a silo with DEM, the location of the front of the rarefaction wave was inferred [194]. Then the formation and propagation of the wave in the discharging process was analyzed. Zaki et al. [195] studied the influence of orifice shape on the discharging process of a flat-bottomed cylindrical silo. Horabik et al. [196] studied the effect of particle shape and filling method on the radial profile of the normal pressure on the bottom of a shallow silo.

3.4.3. Simulation of Grain Separation

Separators were used to separate grains from stones or other impurities. Kannan et al. [197] simulated the separation process of the discarding of a heavy product (stones) and accumulation of a light product (grains) on the deck of a destoner. Then, the effect of deck inclination vibration speed and fluidizing air velocities at the surface on the segregation performance was investigated and their optimal combination value was obtained [198].

Meng et al. [199] conducted simulation research on the separation of whole and broken rice in an indented cylinder separator; the motion trajectory of whole and broken rice with different rotational speed ratios was recorded, and the effect of indent number on the escape angles was investigated. Then the optimal angle of the trough was decided.

4. Research Prospects

DEM could conveniently establish simulation models of agricultural machines and their components, and it was possible to quickly adjust structure and operation parameters within DEM to conduct simulations and collect movement and dynamic data of agricultural materials. Therefore, DEM has become an important approach to help with the design and optimization of the structure and working parameters of agricultural machines and their components. However, the application of research on DEM still remains in its infancy, its application breadth and depth should be developed further.

The DEM modeling accuracy of agricultural materials should be promoted to a higher degree. Most of the existing DEM models of agricultural materials show isotropy and were constituted by spherical particles with the same mechanical properties. However, the construction of agricultural material is usually complicated, and differences exist between different parts of the seed, maize straw cortex and flesh, grape fruit and stem, etc. These differences make agricultural materials isotropic, due to which the physical and mechanical characters remain to be studied further and the models need to be improved accordingly to promote accuracy and reliability, such as the establishment of compound root–soil models. The studies of deformation, failure, and breakage of crop grains and straws are still in a primary stage, and the contact model and model parameters need to be studied further with simulation research of movement and dynamic behaviors, and considering the bending, twisting, cutting, smashing, and damaging behaviors of agricultural materials.

There are also challenges in improving the basic theories of DEM. The contact models in DEM are hypothetical models which approximate the real materials. Certain differences of mechanical relationships exist from the real situation, for example, the adherence between soil particles and tillage equipment still lacks a reliable contact model. To this end, the DEM needs to be improved in its fundamental theories, to make the nonlinear mechanical and dynamic behavior produced by simulations better match the real states. One of the main advantages of DEM is that it is less time consuming, but in the data obtaining process, the mechanical and dynamic information of each DEM particle needs to be calculated. When the DEM particle number is increased, the time spent may increase in a many-fold manner; therefore, how to carry out the calculations more effectively, and improve simulation efficiency is one of the most important research issues.

There is also a need to conduct coupling research along with other methods. DEM is mainly suitable for the research of discrete materials, however, in the interaction process between agricultural materials and agricultural machinery, both discrete materials and continuous mediums should be used. By solid–fluid coupling, complex field situations could be modeled with the consideration of moisture content, air resistance, magnetic force, and so on, with the help of FEM and many-body dynamics. Through solid–solid coupling, the monitoring of the specific resistance of different positions of agricultural machinery components may be realized, and the complex motion of agricultural machinery parts could be modeled; and thus, the agricultural machinery structures could be better designed and optimized.