Discrete element method simulation of disc type furrow openers in paddy soil

For granular materials, discrete element modeling is one of the best computer tools to simulate their behavior and interactions. A field experiment was carried out to evaluate the performance of disc furrow openers in paddy soil. Discrete element simulation was done to develop a 3D DEM model for notched, toothed and double disc furrow openers using EDEM software. Hertz Mindlin contact model with bonding was applied for simulation to fulfill the obligations of the soil moisture and bonding between the cohesive particles. Simulated and field experimental data were compared to determine the applicability in the different working conditions. The results of the simulation validated the applicability of the Hertz-Mindlin contact model with bonding to simulate the no till paddy soil using an extremely narrow tillage tool. The calibrated value of normal and shear stiffness was 5×10 N/m, and the calibrated value of bond normal and shear strength was 3×10 Pa. The relative error (–1.7% to 20.6%) for the double disc furrow opener was lower as compared with that notch typed (29.2% to 44.4%) and toothed type (31.5% to 45.9%) furrow openers.


Introduction 
Physical, structural and biological status of the soil is amended by tillage for better plant growth. Conventional tillage practices degrade the soil and water resources thereby threatening the sustainability of the system [1][2][3][4][5] . In contrast, the conservation tillage system is a low cost production system and resource use efficiency is highest in no tillage system [6] however crop residue management on the soil surface is a difficult task. Numerous types of furrow openers are in use for zero tillage such as disc, hoe and chisel for residue handling and seed placement. Disc type furrow opener is an effective component of seed drill which utilized to develop furrow for seed placement and incorporate the crop residue in the soil layer. Disc opener has an advantage over the hoe type opener due to its low soil disturbance functions [7,8] . The common disc type furrow openers used in the conservation tillage tool are smooth single disc, toothed type, notched type and double disc type. To improve the performance of furrow openers under conservation tillage a number of studies have been executed [9][10][11] . Soil-tool interactions are generally defined in terms of forces developed at the soil-tool boundary and soil particles dislocation [12][13][14] .
Soil-tool interaction modeling has significance to attain the design optimization of tillage tool and operations constraints without performing extensive field tests. Soil tool interaction modeling is relatively complex because of soil variability and discontinues the behavior of soil [15] . Various analytical, empirical and numerical modeling techniques have been applied to study soil tool interaction [16] . Due to the non-homogenous soil structure assumptions of soil failure pattern analytical approach showed low efficiency in simulation various conditions [17][18][19][20][21][22] . To handle the deficits of analytical and empirical methods, numerical modeling approaches have also been used to study the soil tool interaction [15,23,24] . Numerical methods for soil tool interaction can be classified into continuum and discontinuum methods. Continuum numerical methods are further classified as Finite Element Method (FEM), Computational Fluid Dynamics (CFD) and Smooth Particle Hydrodynamics (SPH). Finite element method has been used to study the more complex soil-tool interaction in past research work [24][25][26][27][28] . In FEM simulate the soil as a continuous single object but soil particles move as an independent object during tillage operation thus its assumed that simulation of soil using FEM is limited [28,29] . CFD simulation can be applied for soil tool interactional study (soil dynamic, soil deformation, stress distribution and energy optimization of tool) under various agro-climatic conditions [16,30,31] . However further research is needed to find the model parameters before CFD can be applied to model soil-tool interaction with confidence. SPH is a mesh-free technique, which was developed to resolve the limitations of FEM based methods due to grid distortion [32] . The original SPH method is incompetent to directly resolve the problem of elastic-plastic flows of soil due to the so-called SPH tensile variability [33] . Because of the shortcoming of the analytical method and continuum numerical method, a discontinum numerical method widely known as Discrete Element Method (DEM) is accepted as a potential approach to model mechanical behavior of granular materials. The ability of DEM to focus microstructure level enhances the capacity to understand soil-tool interaction process and leads to improvement of tillage tool design [15] . DEM was first developed for soil mechanics [34] and used to determine the motion between particles. The static contact force between particles and Newton's laws are used to determine the force and motion of particles. The DEM technique can be applied to analyze the mechanical micro and macro-properties of soil.
In the past, most studies were conducted to determine model parameters, evaluation and development of the contact model and tool forces in cohesionless soil using 2D DEM simulation [23,35] . To date, most of DEM simulations were performed for simple earth moving tool (wide blades) in cohesionless soils under bulldozing conditions [36][37][38][39][40] . A few research studies on DEM modeling were conducted for cohesive soil but those also dealt with wide cutting tools [41,42] . Whereas most of the tillage tools have complex shapes and the soils under agricultural production are cohesive with spatial variation in properties. Simulation studies were was conducted using tillage tools such as plow [43] , subsoiler [44,45] and sweep tool [46] in agricultural soils. Mak et al. [47] simulated a simple narrow cutting blade using coarse and fine soil. Most of these studies were conducted using particle flow code (PFC2/3D). But the selection of the appropriate particle contact model, model particle properties and model parameters still a challenge for researchers in DEM modeling.
Proper prediction of soil tool interaction mainly depends on the model parameters [23] . Shmulevich et al. [38] calibrated the model parameters for soil-wide cutting blade using PFC2D. van der Linde [44] applied PFC3D for model parameter determination in sandy soil conditions. Sadek et al. [48] simulated shear behavior of sandy soil and calibrated PFC3D parameters (ball stiffness and bond stiffness) at different moisture content and bulk density level of sandy soil. Coetzee and Els [49] developed and validated a procedure to calibrate the micro material properties required for two dimensional DEM modeling with limitation of granular cohesionless material. Mak et al. [47] found out particle stiffness of coarse and fine soil using PFC3D. Particle contact model also affects significantly DEM simulation results [50] .
The Hertz-Mindline contact model and parallel bond model (PBM) has been recently applied to model the soil tool interaction [40,42,45,47,51,52] . Ucgul et al. [50] made a comparison with the Hysteric Spring Contact Model (HSCM) to Hertz-Mindlin Contact model (HMCM) by simulating wide sweep and determine the model parameters for cohesionless soil using EDEM software. EDEM is multipurpose DEM software with the facility of a user interface [53] . Ucgul et al. [53] used EDEM for 3D DEM tillage simulation to validation of a hysteretic spring contact model for a sweep tool operating in a cohesionless soil. But still it needs to simulate extremely narrow tillage tool like disc type opener in paddy soil using discreet element method especially using EDEM software. The present study was conducted with aimed to develop a discrete element model to simulate disc furrow opener in paddy soil interaction using EDEM.

Field experiment
Field tests were executed at the Experimental Farm, Nanjing Agricultural University, China (32.34°N and 118.36°E) during December 2013 in the post paddy harvested field [14] . Experimental field soil was classified based on international soil textural triangle as clay loam having sand (>0.2 mm) 38.85%, silt (0.002-0.2 mm) 39.84%, clay (<0.002 mm) 21.30% and organic matters 3.18%. Soil plastic limit was 26.7% and the liquid limit was 47.3%. A rice-wheat cropping pattern was followed over a long period of time on the soil of the experimental area. Moisture content and bulk density were determined by the oven drying method of 24 undisturbed soil from the experimental field prior to the commencement tests. Cohesion, internal friction angle, and shear strength were found out using direct shear box apparatus. A digital penetrometer (TJSD-750, Zhejiang Top Instrument Co., Ltd, China) was applied to determine the soil surface penetration resistance (cone index) from ten random locations in the experimental field. At the experiment field, the soil had moisture content 33.3%, bulk density 1.28 g/cm 3 , internal friction angle 12.67° at 0 cm and 8.55° at 10 cm depth, soil cohesion 42.1 kPa at 0 cm and 62.7 kPa at 10 cm depth and average soil cone index 682 kPa.
Figure 1 View of the tested furrow openers Experimental tests were conducted for performance evaluation of three different types of furrow openers (toothed type, notched type and double disc) with external diameter 450 mm in no till paddy soil of furrow openers ( Figure 1). The detail of the key parameters of the furrow opener under investigation is presented in Table 1. A test bench used was equipped with a data acquisition system, power source, tool adjustable frame, tool moving trolley and load cells as detail described in Ahmad et al. [14] . Tillage tools were attached with the depth adjustable frame which had two pairs of parallel connecting strong rods and frame attached to the trolley. As the disc tool moved both the draft force signal data were recorded through the channel data acquisition system using LabVIEW developed model. A completely randomized (3×3×3) factorial experimental design was followed to carry out the study. Three various types of furrow openers (toothed type, notched type and double disc) were tested at three levels working depth (30 mm, 60 mm and 90 mm) and three levels of speed (0.1 m/s, 0.2 m/s and 0.3 m/s). Each test was replicated three times.

DEM simulation 2.2.1 Discrete element model development
Three dimensional models of disc furrow openers with similar dimension to the actual furrow opener used in the field and soil bin with the dimension of 1000 mm×500 mm×200 mm was designed Pro/Engineer (PROE) to use in the DEM simulations. The CAD models were saved in Initial Graphics Exchange Specification (IGES) file type. To minimize the simulation time the size of the soil bin was less as compared to the field.
Reputed computer software for the discrete element simulation, EDEM 2.1 software (DEM-Solutions, Edinburgh UK) having the user interface properties was used for the 3D DEM simulation of the disc furrow openers in paddy soil. The simulation work was executed on a simple desktop computer having CPU 2.2 GHz and 2 GB RAM Capacity, which was significantly lower capacity than the computer (DEL Precision T7500 Intel (R) Xeon (R) CPU X5667 @ 3.07 GHz and 48 GB RAM) used by Ucgul et al. [53] for simulation of sweep tool using similar software. EDEM is a type of user interface software that focuses on basic elements: Particle and geometry (tool, boundary). Particle or a cluster of particles express material particles like soil particles, while boundaries show physical boundaries around the particles, such as a soil bin. Boundaries and other geometric tools can be imported from the CAD files.
An assembly of particles is contained within boundaries in simulating soil-tool interaction.
With the movement of the furrow opener through the particle aggregations, each particle contacts with several other neighboring particles, and the dynamics (displacements and forces) of the particle assembly modifies. EDEM software mainly has three parts as Creator, Simulator and Analyst. The geometry file was imported in the IGES geometry file format.
Because of the bulk of the particles were required to fill the virtual soil bin and to lessening the simulation time in the range of 1-2 h, single phase particle radius 8-10 mm was selected. Hertz Mindlin contact model with bonding was used to fulfill the requirements of the soil moisture and bonding between the cohesive particles. The parameters were selected in the present study according to the constraints of the Hertz-Mindlin contact model with PBM as expressed by EDEM [54] and Ucgul et al. [50] The DEM parameters were classified as material, interaction and bond properties. About similar model parameter values were used by Ahmad et al. [55] The detailed descriptions of these properties are presented in Table 2.

Properties
Value Source

Parallel Bond Configuration
Normal Stiffness/Pa 5×10 7 Mak et al. [47] Shear Stiffness/Pa 5×10 7 Mak et al. [47] Shear strength/Pa 3×10 7 By trial and error Normal strength/Pa 3×10 7 By trial and error 2.2.2 Contact model applied for simulation Several built in models are available in EDEM which are used to develop various interaction between particles to simulate the behavior of various materials, Hertz Mindlin with Bonding contact is one of those which is based in the work of Minddlin [56] combined with the parallel bond contact model (PBCM) which was presented by Potyondy and Cundall [57] are used for the rocky or cohesive materials, in which particles are held together by bonds. van der Linde [44] also pointed out that PBM is relatively much appropriate for agricultural soil simulations. Mak et al. [47] used the PBM contact model for simulation of coarse and fine soil.
Chen et al. [51] also applied the PBM contact model in PFC3D for simulation of sweep in three different types of sandy soils. A cemented like behavior is produced in the bond between adjacent particles. The developed bond can oppose movement (tangential and normal) up to maximum stress (normal and tangential) at which point the bond breaks. Before the bond formation time, particles interact by the Hertz Mindlin contact model [54] .
EDEM's built-in Hertz-Mindlin contact model computes the normal contact forces Fn (N) in a similar manner to Mindlin [56] and Tsuji et al. [58] . 3 (5) where, K hm is the normal stiffness; δ is the normal overlap a; C hm is the normal damping; E is Young's modulus; E e is the effective Young's modulus; R e is the particle radius; V is Poisson's ratio; G is shear's modulus.
The Hertz Mindlin with Bonding contact model can be applied when a finite-sized "glue" bond is required. The tangential and normal motions of the particles in the bond are stopped up to a maximum normal and tangential shear stress. Particles developed bond at formation time of bond. Earlier this moment, the interaction of the particles according to the Hertz-Mindlin contact model. After bonding, the forces (F n , t)/torques (T n , t) on the particle are adjusted to zero and then fixed incrementally every time step according to:  (13) where, A is the area of the bond disk; J is the polar moment of inertia of the disk cross-section; I is the moment of inertia of the disk cross-section about an axis through the contact point; i F, si F are normal and shear component vectors of force; i M, si M are normal and shear component vectors of moment The parallel bond breaks if the maximum shear stress increases over the shear strength or the maximum tensile stress increases over the normal strength [57] .

Selection of the model parameters
Simulation results are affected by the model parameter [23,47,51] . The particle size (ball size) of 8 to 10 mm was used in the simulation to reduce the simulation time. In some studies, the ball size higher than 10 mm was used [47,50,51,53] . Mak et al. [47] selected the ball size 10-20 mm for the determination of the model parameter using PFC3D. Ucgul et al. [53] used the ball size up to 10 mm to determine the appropriate contact model in EDEM software for better discrete element modeling of cohesionless soil. After careful trial and error approach, the bond normal and shear strength was finally assumed 30 MPa. The bond normal stiffness (Kn) and bond shear stiffness (Ks) were selected as 5×10 7 N/m 2 . To achieve the values, several trails were run. In DEM simulation using a narrow tillage tool, the draft force was significantly affected by the bond properties especially the bond shear and normal strength. The thickness of the furrow opener disc (5 mm) was less than the diameter of the particle (10 mm). Thus, the disc mostly moved through the bond area and did not significantly affect the position of the particles. This leads the difference between the simulated and measure draft force of signal disc notched type and toothed type furrow openers.

Comparison of simulated and measured draft force
The DEM simulation was executed using with EDEM using Hertz Mindlin Contact Model with bonding. The 3D views of simulation of notched typed, toothed typed and double disc furrow openers are shown in Figure 2. The 3D simulated view expressed that the particle movement and disturbance was minimum when notched disc furrow opener was operated and this was similar to the field observations. The opener agitates the nearest particle during operation while other particles remained on their position. Figure 2 clearly expressed a comparison of the soil movement with double disc furrow opener quite different from the single notched and toothed type disc type of openers. The particles disturbance area by the double disc furrow opener is significantly higher than the area disturbed by the notched typed and toothed type openers. The cutting width of by double disc furrow opener was 70 mm. This confirms the pervious soil disturbance theory [19] that soil disturbance significantly affected by the width of the tool. During simulation, bond normal stiffness, shear stiffness, bond shear and normal strength significantly affected the simulated results. Chen et al. [51] documented the particles movement of the sandy soil and found that the velocities of the particles are maximum at the soil sweep boundary (from the tip to the shank). It further increased with aggregation of particles on the tool. With increment of time, attached particles adjusted their path behind the sweep and particle velocities of particles ultimately decline to zero. Ucgle et al. [53] determine the soil disturbance area cutting width by sweep in cohesionless soil using DEM simulation. They concluded that 400 mm wide cutting tool developed the 508 mm cutting width and 725 mm thrown widths. The simulation results of furrow openers were compared with the measured results of no till paddy field. The agreements between simulates and measurements were assessed by relative error using following equation: where, R E is relative average error for draft force, %; D stands for draft force; subscripts M and S stand for measurement and simulation respectively. Table 3 shows the relative error for furrow openers.  Figure 3 shows a comparison of measured and simulated results for notch type furrow openers. The simulation results of the draft showed a similar trend of variation with the change of depth and speed to that of measured results. The maximum simulated draft force for notched type furrow opener at 90 mm depth and 0.3 m/s was 497.2 N whereas depth and speed the draft force requirement of notched disc in the paddy field was 762.9 N at 30 mm and speed 0.1 m/s. And at 30 mm of depth with the speed of 0.1 m/s, opener expressed minimum draft force of 134.8 N during simulation which is 29.2% lower than the field results at similar depth and speed. The relative error for the notched opener was between 26.4%-44.4% (Table 3). A similar trend of draft forces in simulation has been observed for the toothed type furrow openers. The toothed type furrow opener also showed the maximum simulated draft force (471.2 N) at 90 mm depth and speed 0.3 m/s (Figure 4). The relative error for the toothed opener was between 20.1%-45.9%.
The comparison of simulated and measured draft force of double disc furrow opener is presented in Figure 5. The lowest relative error (-1.7% to 20.6%) was found for double disc furrow openers as compared to the relative error found in notch type and toothed type furrow openers (Table 3). Simulating the sweep tool in coarse sand using PFC3D, Chen et al. [51] found the relative error between 2% to 31% for soil cutting force at different values (0.50×10 3 -1.75×10 3 ) of particle normal stiffness (Kn). Similarly, Li et al. [45] found the relative error with minimum of 2.96% and maximum 14.95% during the DEM simulation of the subsoiler in agricultural soil. From the results, it clearly shows that the tool having narrow width (single disc, notch type and toothed type) showed higher relative error than that of double disc furrow openers which had relative wider tool width. Therefore, it seems, along with other reasons, the larger particle size (8-10 mm) than the tool width (5 mm) is major reason for higher relative error in notched type and toothed type furrow openers. Increasing the computer power, the size of the particle can be reduced which ultimately can lead to decrease simulation time for each trial. However, the low relative error for double disc furrow openers indicates DEM using EDEM software can simulate the working process of the disc type furrow opener in paddy soil well. The variance in the measure and simulated draft force might also be because of numerous reasons related to simulation and experimental circumstances. In the field, a significant amount of crop residue along with the stubble and roots was present. Therefore, to till the soil under no till environment, the draft force requirement enhanced as the opener had not only cut the soil surface but also chop the residue for its effective penetration into the seed placement zone.
DEM parameters affect the simulation results to large extent [23,47,51] thus the other reason might be due to variation of DEM parameters as because most of the former studies focused wide tillage tool like blades and sweeps under cohesionless soils conditions thus the particle and bond parameters generally available to the wide geometric tool working under non cohesive soil conditions, secondly particle contact model also have a significant effect as Ucgul et al. [50] Ucgul et al. [53] also conducted a study to compare the effect of Hystertic Spring Contact Model (HSCM) and Hertz-Mindlin Contact Model (HMCM) on simulation results and conclude that particle contact model significantly influenced the quality of the results and finally because the capacity of the computer affected the simulation time and thus forced to choose the improper parameter. Therefore, it is required to conduct a study in the future to find out more suitable the particle contact model and their parameters for paddy soil under till and no till condition. The selection and design of appropriate contact model, determination of particle and bond parameters and development of straw and stubble conditions according to the paddy soils are the area where need more comprehensive study for efficient 3D DEM simulation.

Conclusions
Discrete element method was used to simulate the disc openers in paddy soil using EDEM software. The results expressed the ability of Hertz-Mindlin contact model with bonding to simulate the paddy soil using extremely narrow tillage tool. The little higher relative error was due to the parametric effect and larger particle size than the actual soil particle size. However, the partial movement trend was like field tests. With trial and error, these parameters were selected and retrieved the simulated results near to the experimental results with relative error from -1.7% to 45.9%. Therefore, it can be concluded that circular shaped particle with bonding at particle contact points have the ability to produce cohesive soil behavior of agricultural soil however it is required to create aggregative effect during the simulation.
For more improved simulated results, future study will need to consider the relatively similar size of particle as soil and bonding force among the particle due to the moisture in the paddy soil.