Comparative studies of physical models of soil implemented by the method of discrete elements (DEM) on the example of a simple loosening tillage tool

A brief review of the methods used to simulate the interaction of soil-tillage tool is given. The EDEM Academic program was used to practically test the possibilities of using the discrete element method (DEM). A brief description of the physical models of this application used to simulate the interaction of the soil-tillage tool is given. To evaluate the selected physical models (Hertz-Mindlin, Hertz-Mindlin JKR, Hysteretic Spring), a virtual soil channel and a model of a simple tillage tool of rectangular shape have been created. The basic parameters of the soil channel, discrete particles and physical models are described. As a result of the research, two components of the traction resistance vector were obtained, the longitudinal component X and the vertical component Z. Data analysis showed that the average values of forces calculated using physical models of Hysteretic Spring and Hertz-Mindlin JKR differ slightly. With a particle radius of 10 mm, they are equal to X – 597 and 592H; Z – 111 and 97H, respectively. More significant differences when comparing these models with the base model Hertz-Mindlin. For example, in comparison with the physical model of the Hysteretic Spring, the component forces X and Z are less than 18.1%. This is due to the introduction of additional forces imitating the adhesion of particles.


Introduction
To simulate the interaction of a tillage tool with soil, the following numerical methods are used: the finite element method (FEM), computational fluid dynamics (CFD) and the discrete element method (DEM).
The finite element method has found wide application in studies of soil cultivation [1][2][3]. Their analysis shows that the reliability of the forces obtained using this method even when using non-linear soil models is insufficient. The movement of soil aggregates is reproduced either unreliably or very roughly.
Computational fluid dynamics is much less widespread due to significant errors in the forces obtained and the inability to assess the movement of soil aggregates [4][5][6].
The currently most commonly used method for modeling soil cultivation processes is the discrete element method [7][8][9][10]. It can describe soil destruction, deformation, and displacement of soil aggregates. The power characteristics of the process also have high reliability.
Existing studies use a single calibrated physical model. Verification of the results is carried out by comparison with experimental data. At the same time, comparison of data obtained using different physical models is impossible, due to the use of various tillage tools and experimental conditions (diameter and number of particles, speed of movement, size of the soil bin, processing depth). This greatly complicates the process of choosing the most suitable physical model.
For comparative studies, the EDEM Academic program was used, which implements several physical models suitable for modeling soil. This is Hertz-Mindlin, Hertz-Mindlin with JKR Cohesion and Hysteretic Spring.
Hertz-Mindlin (no slip) contact model is the default model due to its accurate and efficient force calculation. In this model the normal force component is based on Hertzian contact theory [11]. The tangential force model is based on Mindlin-Deresiewicz work [12,13]. Both normal and tangential forces have damping components where the damping coefficient is related to the coefficient of restitution as described in [14]. The tangential friction force follows the Coulomb law of friction model as in, for example, [15]. The rolling friction is implemented as the contact independent directional constant torque model [16]. Hysteretic Spring contact model allows plastic deformation behaviors to be included in the contact mechanics equations, resulting in particles behaving in an elastic manner up to a predefined stress. Once this stress is exceeded, the particles behave as though undergoing plastic deformation. The result is that large overlaps are achievable without excessive forces acting upon them, thus representing a compressible material [11]. Hysteretic Spring normal force calculation is based on the Walton-Braun theory described in the following references [17].
Hertz-Mindlin with JKR (Johnson-Kendall-Roberts) Cohesion is a cohesion contact model that accounts for the influence of Van der Waals forces within the contact zone and allows the user to model strongly adhesive systems, such as dry powders or wet materials [11]. In this model, the implementation of normal elastic contact force is based on the JKR theory [18].
Thus, the aim of the work is a comparative study of various physical models describing the interaction between discrete particles and a simple geometry tillage tool under the same experimental conditions.

Physical models
The EDEM Academic software was used for the experiments. The soil was modeled using three different physical models: Hertz-Mindlin, Hertz-Mindlin with JKR Cohesion and Hysteretic Spring. To calibrate the models, we used data obtained in studies of various authors for similar soils [7][8][9][10]. Table 1 shows the global modeling parameters that are used in the Hertz-Mindlin model, as well as the basic parameters in the physical Hertz-Mindlin models with the JKR Cohesion and Hysteretic Spring.  For the Hertz-Mindlin with JKR Cohesion model, it is necessary to additionally set the value of the surface energy (table 3).

Parameters of the virtual soil bin and soil particles
For comparative studies, a virtual soil bin was created, having a length of 2 m and a width of 0.5 m (figure 1). The section of the soil bin used for data collection was chosen as long as 1 m in the central part (figure 2). It was in this area that a steady state was observed according to the results of preliminary experiments. The speed of the tillage tool was taken equal to 1m/s, and the total simulation time was 3s. The data collection interval was 1.02 ... 2.01 s.  The radius of soil particles in various series of experiments was assumed to be 5, 10, and 15 mm. The spherical shape was used exclusively (figure 3). The density of the soil is taken equal to the density of the solid phase of the soil [19]. When backfilling soil particles, the formed soil layer will have a much lower density, which corresponds to the density of natural forest soils. The soil bin was filled using the volume filling function of a static particle factory. The initial height of the soil bin was 0.5 m ( figure 4). After starting the simulation on particle shrinkage stood time interval 0 ... 0.5 s.
The resulting soil layer with a particle radius of 10 mm consists of 31169 particles. The reservoir height is 280 mm, respectively, its volume is equal to 0.28 m 3 and the total mass of particles is 339.5 kg. The density of the soil layer is 1212.5 kg/m 3 , which corresponds to the density of forest soils.

Monitoring of soil dynamic behavior
For visual inspection of the results of modeling and estimation of the displacement of soil particles, color diagrams of particle velocity were used (figure 5). The zero speed of the particles corresponds to blue, the velocity of 1 m/s is red, the intermediate values are shades of green.  The transverse profile of the soil layer was also monitored prior to the passage of the tillage tool (figure 7, a) and after the formation of a furrow ( figure 7, b). For this, in the center of the data removal zone, a transverse soil layer with a thickness of 500 mm was cut off. The resulting profile image was further analyzed in CAD, which allowed us to obtain the exact geometric parameters of the furrow.

Results and discussion
As a result of the study, two components of the vector of traction resistance were obtained -the projection of force on the X axis (longitudinal component) and the projection of force on the Z axis (vertical component) ( figure 8). The average values of forces obtained using physical models of Hysteretic Spring and Hertz-Mindlin JKR differ slightly. With a particle radius of 10 mm, they are equal to X -597 and 592H; Z -111 and 97H, respectively. The difference between the values for X was 0.8%, Z was 12.6%. More significant differences when compared with the base model Hertz-Mindlin. For example, in comparison with the Hysteretic Spring, the component forces X and Z are 18.1% less. This is due to the introduction of additional forces imitating the adhesion of particles.  The change in resistance forces as a function of particle radius was also analyzed. Figure 9, a shows the change in the longitudinal component of the force X. It can be seen that with an increase in the diameter of the particles, the force X increases almost linearly. The accuracy of the linear function approximation for the physical model Hertz-Mindlin is 0.954, Hertz-Mindlin JKR 0.9754 and Hysteretic Spring 0.9876. Figure 10b shows Figure 10 shows sections of the transverse profile of the three furrows. Analysis of their form showed that when using the Hertz-Mindlin basic physical model and the Hysteretic Spring model, the transverse profile of the furrow differs slightly. So, the depth of the furrow relative to the initial surface of the soil layer and the total depth from the bottom of the furrow to the tops of the ridges differ only by 3 mm (47, 50 and 55, 58 mm respectively), and the distance between the centers of the ridges is 17 mm (253 and 236 mm). At the same time, the difference between the physical model Hertz-Mindlin JKR and the basic Hertz-Mindlin are more significant. The depth of the furrow relative to the initial surface of the soil layer differ by 18 mm (50 and 32 mm). This is probably due to an increase in volume due to loosening of the formation. Differences in the depth of the furrow from the bottom to the tops of the ridges are not significant, amounting to 4 mm (58 and 54 mm). The distance between the centers of the ridges differ by 36 mm (236 and 200 mm).  Figure 10. Dependence of the shape of the transverse profile of the furrows on the physical model used.
Analysis of the data showed that the profile that was obtained using the Hertz-Mindlin JKR physical model is the most different, but these differences are not so significant as to draw conclusions about the inadequacy of any of the models for studying the soil profile.

Conclusions
The experiments carried out in the virtual soil bin using the discrete element method showed that data can be obtained in the shortest time possible describing in detail the process of the interaction of the soil -tillage tool system. Quantitative simulation results are significantly affected by both the type of physical interaction and the diameter of soil particles. In this case, the diameter has the greatest effect causing a significant increase in longitudinal and vertical forces. The growth pattern is similar for all models and is close to linear. For example, for the Hertz-Mindlin JKR physical model, with an increase in radius from 5 to 15 mm, component X increased from 330 to 740 N (2.2 times), Y from 52 to 162 N (3.1 times). When modeling it is also possible to study the transverse profile of the formed furrows. The experiment showed that their shape does not have significant differences when using different physical models.
The obtained data are in good agreement on a qualitative level with similar practical and theoretical studies of the interaction of simple tillage tool with the soil [20,21,4,22]. Also, the data obtained using various physical models have good consistency, which confirms the high reliability of the method. However, to obtain more reliable quantitative indicators, it is necessary to conduct additional laboratory studies to calibrate physical models and, based on their results, adjust (if necessary) the input parameters of the physical models.