Simulation of tillage forces and furrow profile during soil-mouldboard plough interaction using discrete element modelling

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Available online at www.sciencedirect.com

ScienceDirect journal homepage: www.elsevier.com/locate/issn/15375110

Research Paper

Simulation of tillage forces and furrow profile during soil-mouldboard plough interaction using discrete element modelling Mustafa Ucgul*, Chris Saunders Agricultural Machinery Research and Design Centre, School of Engineering, University of South Australia, Mawson Lakes, SA, 5095, Australia

article info

The mouldboard plough has a complex shape and its development has mainly relied on

Article history:

resource intensive and time consuming empirical or semi empirical methods. However, if

Received 28 May 2019

the mouldboard plough to soil interaction can be simulated, improved and more energy

Received in revised form

efficient mouldboard ploughs can be designed without performing field tests which may

24 October 2019

only be undertaken at certain times of the year. In this study the interaction between the

Accepted 30 November 2019

soil and a mouldboard plough was simulated using discrete element method (DEM) and results were compared to experimental test and analytical draught force results and furrow profile measurement. It was found that draught and vertical tillage forces can be predicted

Keywords:

within 1.9e16.6% and 0.3e21.5% of the total cutting forces, respectively. This shows that

DEM

DEM can predict draught forces better than analytical method which predicted draught

Mouldboard plough

forces within 3.6e44.4% of the total cutting forces. The comparison of the furrow profile

Draught forces

area results also showed that DEM can predict the furrow profile area within 0.8e16.5% of

Vertical forces

the measured furrow profile area. Results of the study proved that DEM has the potential to

Furrow profile

predict soil-mouldboard plough interaction with a reasonable accuracy.

1.

Introduction

The mouldboard plough is a primary tillage tool and is used to (1) provide soil inversion, (2) create the basis for a seedbed and (3) loosen and aerate the soil. Although the mouldboard plough has some advantages, it is considered as slow and costly due to its perceived inherent high draught force incurred by the way it is used. There are many minimum tillage practices that do not include the use of mouldboard ploughs but provide loosening and creation of seedbeds being

* Corresponding author. E-mail address: [email protected] (M. Ucgul). https://doi.org/10.1016/j.biosystemseng.2019.11.022 1537-5110/© 2019 IAgrE. Published by Elsevier Ltd. All rights reserved.

© 2019 IAgrE. Published by Elsevier Ltd. All rights reserved.

used on farms (Saunders, 2002; Vozka, 2007). However, due to the necessity of inverting soil the mouldboard plough is still a tool of choice, for example, in Australia there is renewed interest in mouldboard ploughing to improve the crop yields of non-wetting sandy soils. Burying the top layer of non-wetting soil and bringing deeper soil to the surface is beneficial for plant growth (GRDC, 2014). Nonetheless, due to (1) increasing energy and labour costs and (2) size of the farms, ploughing should be completed in a short period of time. Thus, higher speeds are generally preferred. According to Ucgul, Saunders and Fielke (2017a) using of higher speeds significantly (1)

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Nomenclature a Ac Afpm Afpp b c dT dw e Ed Ev Efp Fc Fdn Fdt Fdm Fdp Fn Fsn Fst Fvm Ft Fvp Ftm g Hcs He Hl Hmc Hms Hp Hs Ht I kl K1 K2 mr

Indices for sphere or implement Contact area (m2) Measured furrow profile area (m2) Predicted furrow profile area (m2) Indices for sphere or implement Cohesion (kPA) Depth of tip (m) Depth of wing (m) Coefficient of restitution Error between measured and predicted draught force (%) Error between measured and predicted downward vertical force (%) Error between measured and predicted furrow profile area (%) Cohesion force (N) Normal damping force (N) Tangential damping force (N) Measured draught force (N) Predicted draught force (N) Normal total contact force (N) Normal contact force (N) Tangential contact force (N) Measured downward vertical force (N) Tangential total contact force (N) Predicted downward vertical force (N) Total measured force (N) Gravitational acceleration (m s
2) Horizontal cutting side force (kN) Draught force component due to the increase in soil potential energy at the mouldboard (kN) Horizontal landside force (kN) Turning and lifting force (kN) Horizontal mouldboard side force (kN) Horizontal point force (kN) Horizontal share force (kN) Total horizontal force (kN) Moment of inertia (kg m2) Reduction ratio in share width (dimensionless) Stiffness for loading (N m
1) Stiffness for unloading/reloading (N m
1) Rupture distance ratio

reduces the amount of top soil buried into the operation depth and (2) increases the draught force. As a result, design improvement of the mouldboard plough and its components is needed. The mouldboard plough has one of the most complex shapes of all tillage implements and its development has been mainly due to empirical or semi empirical methods (e.g. Arvidsson & Keller, 2011; Desbiolles, Godwin, Kilgour &Blackmore, 1997; Li, Lobb & Lindstrom., 2007; Mari et al., 2015; Saunders, 2002). Although empirical methods provide practical information; performing experiment for each field (e.g. soil type, low or high cohesive soil) and operation condition (e.g. speed, depth, share width) is time

m meq M Mr nc nk Ng Nc Na r rcon req R Uabn Uabt U0  Uabn  Uabt U¨ Y Vn w wT ww

59

Mass (kg) Equivalent mass (kg) Moment (N m) Moment due to rolling friction (N m) Damping factor Stiffness factor Dimensionless number for soil density Dimensionless number for cohesion Dimensionless number for adhesion Radius (m) Perpendicular distance of contact point from the centre of mass (m) Equivalent radius (m) Rotational acceleration (rad s
2) Normal component of the relative displacement (m) Tangential component of the relative displacement (m) Residual overlap (m) Normal component of the relative velocity (m s
1) Tangential component of the relative velocity (m s
1) Translational acceleration (m s
2) Yield strength (MPa) Velocity (ms
1) Width of mouldboard plough (m) Width of tip (m) Width of wing (m)

Greek letters a Rake angle ( ) Tip rake angle ( ) aT Share rake angle ( ) aw b Share approach angle ( ) d Angle of soil to metal friction ( ) g Bulk unit weight of soil (kN m
3) 4 Internal friction coefficient of soil ( ) q Mouldboard angle ( ) m Coefficient of friction Coefficient of rolling friction mr Unit vector of angular velocity lq x Cohesion energy density (J m
3)

consuming and costly. In addition, manufacturing new mouldboard designs is also expensive. If the mouldboard plough to soil interaction can be modelled using computers, costly and time-consuming manufacturing and field-testing stages can be reduced. Analytical methods (e.g. Saunders, Godwin & O’Dogherty, 2000; Godwin, O’Dogherty, Saunders & Balafoutis, 2007) and finite element methods (FEM) (e.g. Bentaher et al., 2013; Formato, Faugno & Paolillo, 2005; Zhang, Cai, Wang, Zhang, & Liu, 2018) are commonly employed to model soilmouldboard plough interaction. Due to their quasi-static or dynamic condition assumptions; analytical models only

60

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

examine soil failure forces and cannot model the movement of the soil. Although there are some advantages in using FEM, the assumption of continuity is not always valid, hence a change in the soil structure or soil translocation cannot be predicted (Asaf, Rubinstein, & Shmulevich, 2007). A method that overcomes the shortcomings of empirical/ semi-empirical, analytical and finite element methods is discrete element method (DEM). DEM is a discrete numerical method which can be used to model soil-tool interaction (e.g. Asaf, Rubinstein & Shmulevich, 2006; Ucgul, Saunders & Fielke, 2017b). DEM is based on modelling the individual contacts (physical interactions) between particles where many particles are used to represent a soil bulk. As soil media consists of many small particles, the main limitation of DEM when modelling soil is the number of contacts between the particles hence enormous computation time. In the past computer technology was not suitable for this approach to considered, but with the recent rapid development of software and hardware technology DEM is now being used to model granular materials. However, as actual soil particles are too small (<2 mm), using of real particle sizes is currently still not viable. Therefore, larger than actual particle sizes are generally used in the simulations. Using larger than actual particles sizes has been employed in several studies (e.g. Mak, Chen & Sadek, 2012; Barr, Ucgul, Desbiolles & Fielke, 2018;  rez, Cueto, & Ramon, 2014). In DEM a Bravo, Tijskens, Sua suitable contact model and an accurate calibration process is required to obtain an accurate simulation. Reliability of the DEM results are related to the accuracy of the contact model (Ucgul, Fielke & Saunders, 2014a; Yan, Wilkinson, Stitt & Marigo, 2015). Basic elastic contact models generally provide less accurate results while more complex plastic models can generally provide more accurate results (Ucgul et al., 2014a). A study carried out by Ucgul, Fielke, and Saunders (2015) showed that a linear cohesion integrated hysteretic spring contact model can be used to accurately model soil-tool interaction. In the model, contact forces between the particles are calculated using a hysteretic spring model while cohesion is defined by adding a constant cohesive force between the particles. The cohesive force is applied through the centre of the DEM particle. DEM modelling of soil to a oneethird scale mouldboard plough interaction has been carried out by Ucgul et al. (2017b). The interaction between soil and a full-scale mouldboard plough was also modelled by Ucgul et al. (2017a) to investigate the top soil burial. However, no DEM study has been conducted to model both tillage forces and furrow profile for soil-mouldboard plough interaction in terms of different speeds, depths and mouldboard share widths. In the experimental work of Saunders (2002) the effect of speed, operation depth and mouldboard share width on the tillage forces and furrow profile were investigated. An analytical model was also developed to predict the draught forces. In the current study the work of Saunders (2002) was simulated using DEM and the force and furrow profile results were compared to experimental results. The DEM predicted draught force results were also compared to analytical results of Saunders (2002). The specific aim was to develop a calibrated DEM model to examine the effects of width, speed and depth on tillage forces and furrow profiles.

2.

Methodology

2.1.

Soil bin tests

To investigate the interaction between soil and the mouldboard plough a series of experiments were performed by Saunders (2002). The experiments were performed under the controlled laboratory conditions using an indoor soil bin facility at Cranfield University, Silsoe, UK. A sandy loam soil (68.1% sand, 22.1% silt and 9.8% clay) was used in the experiments with an average moisture content and bulk density of 10.85% (Coefficient of variation (CV) ¼ 2.8%) and 1530.1 kg m
3 (Coefficient of variation (CV) ¼ 1.7%), respectively (both dry basis). Internal friction coefficient of the soil was measured as 30.75 whereas the soil metal friction was measured as 24 . Internal friction coefficient of the soil was measured using triaxial test whereas soil metal friction was measured using direct shear test. The cohesion was also measured as 10,540 Pa. The adhesion of the soil bin soil was negligible. A commercial Kverneland™ No 8 mouldboard assembly, including point, share and landslide was used in the experiment to investigate draught and vertical forces along with soil profile at different operation widths, speeds and depths. The mouldboard plough was fitted to a two-dimensional Extended Octagonal Ring Transducer (EORT), which was linked to signal processing and a data logger recorded the measured values at 103 Hz, to enable the draught and vertical forces to be measured. Calibration was carried out on the EORT using known static loads to match the gauges to the signalprocessing unit. Experiment parameters which were used in the soil bin tests are given Table 1. Each test was replicated three times and the experiment carried out in a totally randomized order. Due to the nature of the mouldboard plough action only one run was possible in each bin preparation. To simulate a true situation an open furrow was produced at one side of the soil bin, which gave the plough assembly an open furrow into which to turn the soil slice. The false furrow was removed to the same depth as the depth of test that was to follow. Into this false furrow a pre-furrow was ploughed using the mouldboard at the corresponding test depth to provide the same open furrow profile as that left by the plough in a field situation. This was needed to ensure that the soil slice being turned had an open furrow to fall into, as in a field situation.

Table 1 e Experiment parameters (Saunders, 2002). Test no 1 2 3 4 5 6 7 8 9 10 11

Width (mm)

Depth (m)

Speed (ms
1)

355 457 508 508 508 508 508 508 508 508 508

0.175 0.175 0.125 0.175 0.225 0.125 0.175 0.225 0.125 0.175 0.225

1.250 1.250 1.250 1.250 1.250 2.083 2.083 2.083 2.778 2.778 2.778

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

This would ensure no spurious forces were encountered by the soil being pushed into soil in the previous furrow or the side of the soil bin (Saunders, 2002). After the pre-furrow was generated the mouldboard plough was run and the measurements were taken. The furrow profile was measured using a pin type 1 m width profile meter (Saunders, 2002) (Fig. 1) which gave an indication of the distance the soil was moved as a result of each plough setting.

61

In the model turning and lifting force (Hmc) was given as; Hmc ¼ (g/g),(wT dT þ wwdw) , Vn , (1-(1-(sinq tand)cosq)

(5)

Force due to potential energy (He) was also expressed as; He ¼ g ,(wT dT þ wwdw) , (mrdw)

(6)

Landside drag force was calculated as;

2.2.

Analytical force prediction Hl ¼ (Hcs þ Hms) tand

To predict the draught forces acting on a mouldboard plough an analytical method was developed Saunders (2002). In the method, forces acting on the mouldboard plough were divided into horizontal point force (Hp), horizontal share force (Hs), turning and lifting force (Hmc), force due to potential energy (He) and horizontal landside force (Hl) (Fig. 2). Horizontal point force (Hp) was calculated as; Hp ¼ [[(gdT2Ng þ cdTNc),(wT þ 0.5dT(mr-(1/3(mr-1)))] þ [((w2NadT)/g),(wT þ 0.6dT)]],sin(aT þ d)

(1)

(2)

Hcs ¼ [[[(gdw2Ngw þ cdwNc) þ ((w2Nadw)/g)],(ww-kldw)],sin(aw þ d)],(cosb)tand (8)

Hms ¼ (g/g),(wT dT þ wwdw) , Vn ,sinq(1-sinq tand) tand

Na ¼ (tanb þ cot (b þ 4)) / ((cos (a þ d) þ sin (a þ d) cot (b þ 4)) , (1 þ tanb cota)) (3) where b was defined as; (4)

where a is the rake angle of the mouldboard plough and is defined as the angle between the implement’s plane and the horizontal plane in the direction of travel.

Fig. 1 e Profile meter used in the experiments by Saunders (2002).

(9)

The total draught force (Ht) is then calculated as; Ht ¼ Hp þ Hs þ Hmc þ He þHl

(10)

Details of the analytical model can be found in Saunders (2002).

2.3.

where Na in both Eqs. (1) and (2) was defined as;

b ¼ arctan (1 / (mr ecot a))

where Hcs and Hms are horizontal cutting side force and horizontal mouldboard side forces, respectively. Hcs was defined as;

while Hms was expressed as;

where mr is the rupture distance ratio and is defined as the forward rupture distance divided by the depth of crescent failure. Horizontal share force (Hs) was defined as; Hs ¼ [[(gdw 2Ng þ cdwNc) þ ((w2Nadw)/g))],(ww þ kldw)], sin(aw þ d)

(7)

Discrete element method (DEM) simulation

The DEM simulation was carried out using EDEM 2.7™ software operating on a DELL Precision T7610 Dual Xeon E5-2680 v3 @ 2.5 GHz and 128 GB RAM computer. EDEM 2.7™ software was developed by DEM Solutions Ltd. in Edinburgh, Scotland (UK). In order to simulate soil mouldboard plough interaction following steps were undertaken: (a) the material and interaction properties of the particles were entered in to the software, (b) the CAD models of the geometries (i.e. mouldboard plough and the soil-bin) were imported, (c) a particle factory, in where the particles are generated (EDEM, 2011), was created (to generate the particles for the simulation) and (d) simulations were carried out. EDEM 2.7™ provides accurate solutions. It also has extensive post-processing tools for creating graphs, videos and JPG images of the simulation which help users to analyse the simulation in detail. By using EDEM 2.7™ the soil bin experiments were replicated. A linear cohesion integrated hysteretic spring contact model suggested by Ucgul et al. (2015) was used to model the soil-mouldboard interactions. The model allowed the particles to behave in a linear elastic manner up to a predefined stress and when the total stress on the contact area exceeded the predefined stress (which is the yield strength) in the model, the particles behaved as though undergoing plastic deformation. The cohesion between the particles was defined by directly adding a cohesion force to the normal contact forces. In the hysteretic spring contact model the total normal (Fn) and tangential (Ft) forces were calculated as;

62

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Fig. 2 e Diagram of the components of the draught force acting on the plough (Godwin et al., 2007).

Fn ¼ Fns þ Fnd

(11)

Ft ¼ Fts þ Ftd

(12)

where Fns, Fnd, Fts and Ftd are the normal contact force, normal damping force, tangential contact force and tangential damping force, respectively. Normal contact force was calculated as per EDEM (2011) as; 8 < K1 $Uabn Fns ¼
K2 $ðUabn
U0 Þ : 0

loading unloading=loading unloading

(13)

where Uabn is the normal component of the relative displacement, Uo is the residual overlap. K1 and K2 are the loading and the unloading stiffnesses, respectively. As per Walton (2006), K1 was calculated as; K1 ¼ 5 req min(Ya, Yb)

(14)

where Y is the yield strength and req is the equivalent radius and defined as EDEM (2011), 1/req ¼ 1/ra þ 1/rb

(15)

where r is the radius for the individual particles a and b. Following Walton and Braun (1986), K2 was calculated as; K2 ¼ K1/e

2

(16)

where e is the coefficient of restitution. The residual overlap was updated in each time step as; 8 < Uabn$ð1
ðK1 =K2 ÞÞ U0 ¼ U0 : Uabn

loading unloading=loading unloading

The tangential contact force was calculated as per EDEM (2011) as; Ft s ¼ -nk K1 Uabt

(18)

where Uabt is the tangential component of the relative displacement. nk is the stiffness factor which was defined as the ratio of tangential stiffness to normal loading stiffness. The normal and the tangential damping forces were calculated using; Uabn Fnd ¼ -nc ((4 meq K1)/(1 þ (p/ ln e)2))
1/2 

(19)

Ftd¼- ((4 meq nk K1)/(1 þ (p/ ln e)2))
1/2  Uabt

(20)

Uabt are the normal and tangential compowhere  Uabn and  nents of the relative velocity, respectively. nc is the damping factor which controls the amount of velocity dependent damping. meq is the equivalent mass and defined in EDEM (2011) as; 1/meq ¼ 1/ma þ 1/mb

(21)

where m is the mass for the individual particles a and b. The total tangential force was limited to the lessor of either the calculated tangential force or the sliding friction force; Ft ¼ -min (nk K1 Uabt þ Ftd, m Fns)

(22)

The magnitude of the moments caused by total tangential force (M) and the rolling resistance (Mr) were calculated following Raji (1999) by: M ¼ rcon Ft

(23)

Mr ¼ -mr Fn s rcon lq

(24)

(17)

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

63

where rcon is the perpendicular distance of the contact point from the centre of mass, mr is the coefficient of rolling friction and lq, is the unit vector of angular velocity at the contact point. After interacting with other particles, the new position of a particle was calculated by integrating Eqs. (25) and (26). U¨ ¼ (Fn þ Ft)/m

(25)

R ¼ (M þ Mr)/I

(26)

To model soil cohesion, a cohesive force was added to the total normal forces. The inter-particle friction was assumed to restrict the tangential element motion in the governing equations, thus the cohesion force was not added in the tangential direction. The magnitude of the cohesion force was calculated as (EDEM, 2011); Fc ¼ x Ac

(27)

where x is the cohesion energy density which is defined as the energy needed to remove a particle from its nearest neighbours divided by the total volume of the removed particle and Ac is the contact area. After adding the cohesion force, Eq. (11) becomes; Fn ¼ Fns þ Fnd þ Fc

(28)

A 20 m long x 2 m wide x 0.6 m deep virtual soil bin was created using 20 mm diameter spherical particles. No particle size distribution has been given in Saunders (2002). So, particles were randomly generated in the range of 0.75e1.5 times the nominal size (between 15 and 30 mm diameter). Similar particle sizes have been used in the literature for simulating soil-tool interaction (i.e. between 10 and 20 mm by Mak et al. (2012); between 19 and 21 mm by Ucgul et al. (2017a)). This particle size was chosen based on the available computational power. Particle sizes smaller than this caused long computational times. Data obtained from the simulation results has shown that the energy dissipated with rotational movement accounts for 0.021% of total energy. This shows that spherical particles are appropriate for the simulations. The bulk density of the test soil (1530 k gm
3) was achieved using the following method. Firstly, 33,048 kg DEM particles (in total 2,662,670 particles) were generated in the virtual soil bin. Then the generated particles were pushed down using a virtual plane by 60 mm. This gave the bulk density of the test soil (33,048 kg/(20 m  2 m  0.54 m) ¼ 1530 kg m
3). Then the CAD model of the mouldboard plough (for each 355, 457 and 508 mm width) was imported into EDEM 2.7 ™ and mouldboard speeds and depths were set to represent the depths and speeds used in the experiments (as given in Table 1). Subsequently a false furrow was removed to the same depth as the depth of test and then a pre-furrow was generated. After generating the pre-furrow the tool was run again and force and soil profile measurements were taken. The draught and vertical forces were taken over the mid-section of the soil bin where they were stable (Fig. 3a). The furrow profile was also measured in this stable section. After cutting a slice in the soil bin, the crossesection profile of the furrow was obtained (Fig. 3b). Thereafter the particles of the boundary particles

Fig. 3 e DEM simulation (a) Force measurement (b) Simulated furrow profile. were selected and their coordinates were exported to Microsoft Excel for generating furrow profiles. After the soil profiles were generated, the profile curves obtained from soil bin tests, by Saunders (2002) and those predicted, were fitted. The measured and predicted curves were superposed in Microsoft Excel based on from their starting points and depths as the starting point and ploughing depth of each profile are known. As larger than actual particle sizes are used to model soil media, the DEM parameters must be calibrated (so that these larger particles can mimic the actual soil bulk behaviour). The DEM parameters used in the study are shown in Table 2. Some of the parameters were provided by Saunders (2002) while others were taken from the literature or calibrated. In Table 2 cohesive energy density, which is defined as the energy needed to remove a particle from its near neighbours divided by the total volume of the removed particle (J m
3 ¼ N m
2) (EDEM, 2011), was assumed to be equal to the cohesive strength measured by Saunders (2002). Due to lack of test soil, the calibration process was based on matching simulation results to actual measured results of draught and vertical forces at 0.225 m operation depth and 2.778 m s
1 tool speed for 508 mm width mouldboard plough (which were measured by Saunders (2002)). This method was also used by Chen, Munkholm, and Nyord (2013). Using the known parameters, either provided by Saunders or available in the literature (given in Table 2), and varying the coefficient of restitution between soil and soil and yield strength of the soil, smallest relative error to the measured total force (Eq. (31)) was targeted using a trial and error approach (Table 3). After the

64

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Table 2 e DEM parameters used in the simulations. Property

Value

Density of sand particles (kg m
3) Density of steel (kg m
3) Shear modulus of soil (Pa) Shear modulus of steel (Pa) Poisson’s ratio of soil Poisson’s ratio of steel Yield strength of the soil (Pa) Coefficient of restitution of soilesoil Coefficient of friction of soilesoil Coefficient of friction of sand-steel Cohesive energy density between soilesoil (J m
3) Adhesion between soil-tool (Pa)

Source

2600

Huser and Kvernvold (1998)

7861 5  107

Hudson Tool Steel (2016) Academia (2015)

7.9  10

10

Hudson Tool Steel (2016)

0.3 0.3 2  106

Asaf et al. (2007) Budynas and Nisbett (2012) Calibrated

0.3

Calibrated

0.6

Saunders (2002)

0.45

Saunders (2002)

10,540

Saunders (2002)

0

Saunders (2002)

calibration the closest matching draught and downward vertical forces were 5.664 kN and
0.802 kN. It can be said from the calibration results that draught and vertical downward forces were calibrated with a reasonable accuracy.

3.

Results and discussion

3.1.

Tillage forces

A comparison of soil bin test, DEM simulation and analytical force prediction results for different width mouldboard ploughs at 0.175 m operation depth and 1.250 ms
1 speed are given in Fig. 4. The measured values are the mean values reported by Saunders (2002). The coefficient of variation of the measured draught forces were 5.7%,14.6%, 3.6% for 355, 407 and 508 mm widths, respectively. A comparison of soil bin test, DEM simulation and analytical force prediction results for 508 mm width mouldboard plough in terms of three different depths (0.125, 0.175 and 0.225 m) and speeds (1.250, 2.083 and 2.778 ms
1) is also given in Fig. 5 (a), (b) and (c). The coefficient of variation of the measured draught forces were between 0.6% and 14.6%.

Table 3 e Results of calibration procedure. Coefficient of restitution 0.1 0.3 0.5 0.7 0.3 0.3 0.3 0.3

Yield strength (Pa) 2 2 2 2 5 1 2 4

106 106 106 106 105 106 106 106

Simulated draught force (kN)

Simulated vertical downward force (kN)

Simulated total force (kN)

Relative error (%)

4.351 5.664 5.378 5.056 5.132 5.304 5.664 5.768


0.799
0.802
0.719
0.811
0.965
0.865
0.802
0.700

4.425 5.721 5.426 5.121 5.222 5.374 5.721 5.810


21.1 2.0
3.3
8.7
6.9
4.2 2.0 3.6

Measured total force is 5,609 kN.

Fig. 4 e Measured, DEM simulated and analytically predicted tillage forces for different width mouldboard ploughs (355, 407 and 508 mm) at 1.250 ms¡1 speed and 0.175 m operation depth.

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

65

Fig. 5 e Measured, DEM simulated and analytically predicted tillage forces for 508 mm width mouldboard plough at three different speeds for (a) 0.125 m (b) 0.175 m (c) 0.225 m operation depths.

Figure 4 shows that draught force increases with the increase of the furrow cutting width while increasing of the tool cutting width does not have any considerable effect on downward vertical force. Similarly, Fig. 5 shows that draught forces increases with the increase of speed and depth whereas speed does not have any major effect on downward vertical force. It was also found that increasing

ploughing depth increased draught forces more than ploughing speed. Test results showed that increasing of ploughing speed from 1.250 m s
1 to 2.778 m s
1 increased draught force between 3 and 21% while increasing of ploughing depth from 0.125 to 0.225 m increased draught force between 56 and 77%. DEM simulation results also confirm the test results.

66

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Results of this study are also matching with the test results performed by Ibrahmi, Bentaher, Hamza, Maalej, and Mouazen (2015). When it is compared to the FEM analysis carried out by Ibrahmi et al. (2015) who used the linear form of the extended DruckerePrager model, it can be said that DEM predicts draught forces better than FEM in particularly for increasing speed. When the DEM simulation results were compared to analytically predicted results of Saunders (2002) it was found that more accurate draught force prediction can be achieved using DEM. Error between the measured and predicted draught (Ed) and downward vertical (Ev) forces were evaluated in terms of total measured force (Ftm) as; Ed ¼ |Fdm-Fdp| Ftm
1

(29)

Ev ¼ | Fvm-Fvp| Ftm
1

(30)

where Fdm and Fdp are measured and predicted draught (DEM or analytical) forces, respectively. Fvm and Fvp are measured and predicted downward vertical forces, respectively. Ftm is the total measured force and was calculated as; Ftm¼ ((Fdm)2 þ (Fvm)2)0.5

(31)

Coefficient of variation values of the measured downward vertical forces were high (between 16.5% and 55.8%). In addition, the magnitudes of the downward vertical force were smaller than draught force values therefore using standard relative error method is not suitable to evaluate the error prediction. So, error between the measured and predicted draught (Ed) and downward vertical (Ev) forces were evaluated in terms of total measured force (Ftm). Similar methods were also used in the literature (i.e. Chen at al, 2013; Ucgul et al., 2017b). The calculated errors between measured and predicted draught (for both DEM and analytical prediction) and vertical forces are presented in Table 4. The regression plots of measured vs DEM predicted and measured vs analytically predicted draught forces are also shown in Fig. 6a and b. As it can be seen from Table 4 and Fig. 6, a smaller relative error and better correlation were obtained for draught force prediction when DEM is used. The deviation between measured and DEM simulated tillage forces can be attributed to the (1) larger than actual particles sizes used in the simulations, (2) limits of calibration for DEM parameters, and (3) errors during the experimental process (e.g. inconsistencies between the soil bin preparation).

Table 4 e Errors between measured and predicted draught force (Ed) and measured and predicted vertical downward force (Ev). Test no

Width (mm)

Depth (m)

Speed (ms
1)

Ed (%) DEM

Ed (%) Analytical

Ev (%) DEM

355 457 508 508 508 508 508 508 508 508 508

0.175 0.175 0.125 0.175 0.225 0.125 0.175 0.225 0.125 0.175 0.225

1.250 1.250 1.250 1.250 1.250 2.083 2.083 2.083 2.778 2.778 2.778

5.80 7.00 11.6 11.7 7.70 16.6 4.10 2.60 6.60 1.90
3.00

36.3 27.6 44.3 33.3 22.9 42.4 22.8 14.6 29.9 13.9 3.60


13.9 4.00 0.30
8.40
8.90
1.00
11.6
12.2
16.8
21.5
5.50

1 2 3 4 5 6 7 8 9 10 11a a

Calibrated dataset.

Fig. 6 e Correlation between measured draught force vs (a) DEM predicted draught force and (b) Analytically predicted draught force.

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

67

Fig. 7 e Measured and DEM predicted soil profiles for (a) 355, (b) 405 and (c) 508 mm widths at 1.250 ms¡1 speed and 0.175 m depth and (d) 508 mm widths at 2.778 ms¡1 speed and 0.175 m depth.

68

3.2.

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Furrow profile

The measurement of the soil profiles was an important part of the experiment; this indicates the soil movement from each of the plough settings. The soil profile was randomly selected in the steady state region in where draught forces are stable as it was performed by Saunders (2002). Measured and DEM simulated soil profile measurements were presented in Fig. 7. Figure 7a, b and c show that there is no difference in the furrow profile for the three widths of operation where speed is maintained at 1.250 m s
1 and depth is maintained at 0.175 m in both experiment and DEM simulation. It was also found that increasing the share cutting width did not increase the width of the furrow bottom (which was between 0.4 and 0.5 m). Conversely, a very small increase (about 0.1 m) in the furrow bottom width from the slowest speed to the fastest was observed at each operation depth (Fig 7c and d). It was also observed that the profiles have moved further from the furrow wall at the top of the profile where there is less obstruction from the previous furrow. At the slower speeds, the dip between the furrow profile measured and the previous furrow produced in the soil bin is noticeable, where the profile peaks between 0.5 and 0.7 m from the furrow wall, depending on depth. As the speed is increased, the dip between furrows tends to disappear and the profiles have levelled off after one metre (Fig 7c and d). As seen in Fig. 7, DEM predicted open furrow area between 0 and 0.1 m is less than that of the test results. This can be attributed to the larger than actual particle sizes used in the simulations. As the cutting-edge thickness is much smaller than the particle size, instead of cutting the plough blade pushes particles and this caused particles to fall into the furrow profile due to the gravity overcoming cohesion and friction. The quantitative comparison between the measured (Afpm) and DEM simulated (Afpp) furrow areas were carried out by considering the furrow area under the “x” axes. The area under the x axis was calculated using the trapezoid rule (numerical interaction). If we consider the coordinates of two data points as (x1, y1) and (x2, y2), the area between two data points were calculated as ((y1þy2)/2)*(x2-x1). The total furrow area below the curve was calculated as the total area of all trapezoids (areas between two data points) below the curve. Error between the measured and predicted furrow profiles (Efp) was evaluated as; Efp ¼ |Afpm-Afpp| / Afpm

(32)

where Afpm and Afpp are measured and predicted furrow profiles, respectively. The calculated errors between measured and predicted furrow profiles are given in Table 5. Results showed increasing ploughing depth increased furrow area more than ploughing speed. The correlation between the measured and DEM predicted furrow profiles was also given in Fig. 8. As shown in Table 5 and Fig. 8 correlation between measured and predicted furrow areas are quite high. Note that the coefficient of variation for the DEM predicted furrow area of calibrated dataset was calculated as 5.9%. This shows that the variability of the furrow profile is quite low. Similar test results were also observed by Ibrahmi et al. (2015) in the case of ploughing depth. While a more remarkable increase in furrow area was measured by Ibrahmi et al.

Table 5 e Error between measured and simulated furrow profiles (Efp). Test no 1 2 3 4 5 6 7 8 9 10 11a a

Width (mm) 355 457 508 508 508 508 508 508 508 508 508

Depth Speed (m (m) s
1) 0.175 0.175 0.125 0.175 0.225 0.125 0.175 0.225 0.125 0.175 0.225

1.250 1.250 1.250 1.250 1.250 2.083 2.083 2.083 2.778 2.778 2.778

Afpm (mm2)

Afpp (mm2)

Efp (%)

61,253 59,002 36,932 57,125 69,079 45,737 62,179 71,482 46,464 66,362 76,487

66,928 60,894 42,579 56,721 80,438 45,387 60,799 81,873 46,201 63,688 81,915


9.3
3.2
15.3 0.7
16.4 0.8 2.2
14.5 0.5 4.0
7.1

Calibrated dataset.

Fig. 8 e Correlation between measured furrow area vs DEM furrow area.

(2015) in the case of speed. This might be attributed to the method of measuring the furrow area by Ibrahmi et al. (2015) who measured the furrow area after removing the loose soil in the furrow area. The difference between measured and DEM simulated furrow profiles can be attributed to the (1) larger than actual particle sizes used in the simulations (using smaller particle sizes was suggested by Ucgul, Fielke, and Saunders (2014b) and Wang et al. (2019)), (2) limits of calibration for DEM parameters, (3) errors during the experimental process (due to profile meter) and (4) profile variability along the furrow.

4.

Conclusion

In this study, a calibrated DEM model was developed to simulate the tillage forces and furrow profiles measured by Saunders (2002) for soil-mouldboard plough interaction. Predicted draught forces were also compared to the analytical draught force prediction of Saunders (2002). Both test and DEM

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

predicted results have shown that (1) increasing the width of mouldboard plough causes only minimal increase in the width of the furrow profile, (2) increasing ploughing depth increased draught forces more than ploughing speed and (3) increasing ploughing speed reduces the height of the furrow profile and fills the dip between the furrow profile and the previous furrow profile. The study has also proven that the developed DEM model simulated well the measured draught and vertical forces and furrow profiles. And when it is compared to the analytical method, better draught force prediction can be achieved using DEM. This study shows that when its parameters are calibrated, DEM overcomes the shortcomings of empirical/semi-empirical, analytical and finite element methods and can simulate the tillage forces and furrow profiles. To improve the results future work will need to consider (1) developing an improved contact model, (2) developing an improved calibration procedure with using a statistical approach (such as design of experiment) which considers the effects of all micromechanical parameters in the inter connections and (3) using closer to actual particle sizes. The current work will also be extended to simulate the pressure distribution over the mouldboard surface.

Acknowledgments The research has been funded by Australia’s Grains Research and Development Corporation (GRDC) and Department of Primary Industries and Regional Development (DPIRD) DAW 00244. The authors would also like to acknowledge the staff and facilities at Cranfield University, Silsoe, UK for support of the experimental soil bin testing.

references

Academia. (2015). Some useful numbers for rocks and soils. Retrieved from http://www.academia.edu/4056287/Some_Useful_ Numbers_for_rocks_and_soils. Arvidsson, J., & Keller, T. (2011). Comparing penetrometer and shear vane measurements with measured and predicted mouldboard plough draught in a range of Swedish soils. Soil and Tillage Research, 111(2), 219e223. Asaf, Z., Rubinstein, D., & Shmulevich, I. (2006). Evaluation of link-track performances using DEM. Journal of Terramechanics, 43(2), 141e161. Asaf, Z., Rubinstein, D., & Shmulevich, I. (2007). Determination of discrete element model parameters required for soil tillage. Soil and Tillage Research, 92(1e2), 227e242. Barr, J. B., Ucgul, M., Desbiolles, J. M., & Fielke, J. M. (2018). Simulating the effect of rake angle on narrow opener performance with the discrete element method. Biosystems Engineering, 171, 1e15. Bentaher, H., Ibrahmi, A., Hamza, E., Hbaieb, M., Kantchev, G., Maalej, A., et al. (2013). Finite element simulation of moldboardesoil interaction. Soil and Tillage Research, 134, 11e16.  rez, M. H., Cueto, O. G., & Ramon, H. Bravo, E. L., Tijskens, E., Sua (2014). Prediction model for non-inversion soil tillage implemented on discrete element method. Computers and Electronics in Agriculture, 106, 120e127.

69

Budynas, R. G., & Nisbett, K. J. (2012). Shigley’s mechanical engineering design. McGraw-Hill Education. Chen, Y., Munkholm, L. J., & Nyord, T. (2013). A discrete element model for soilesweep interaction in three different soils. Soil and Tillage Research, 126, 34e41. Desbiolles, J. M. A., Godwin, R. J., Kilgour, J., & Blackmore, B. S. (1997). A novel approach to the prediction of tillage tool draught using a standard tine. Journal of Agricultural Engineering Research, 66(4), 295e309. EDEM. (2011). EDEM theory reference guide. Edinburgh, UKCopyright©: DEM Solutions Ltd. Formato, A., Faugno, S., & Paolillo, G. (2005). Numerical simulation of soil-plough souldboard interaction. Biosystems Engineering, 92(3), 309e316. Godwin, R. J., O’Dogherty, M. J., Saunders, C., & Balafoutis, A. T. (2007). A force prediction model for mouldboard ploughs incorporating the effects of soil characteristic properties, plough geometric factors and ploughing speed. Biosystems Engineering, 97(1), 117e129. GRDC (Grain Research and Development Corporation). (2014). Soil inversion (mouldboard ploughing). Retrieved from: http://www. grdc.com.au/Research-and-Development/Major-Initiatives/ Non-wetting-soils/Soil-inversion-mouldboard-ploughing. Hudson Tool Steel. (2016). P20 Mold steel. Retrieved from http:// www.hudsontoolsteel.com/technical-data/steelP0. Huser, A., & Kvernvold, O. (1998). Prediction of sand erosion in process and pipe components. In BHR group conference series publication (Vol. 31, pp. 217e228). Mechanical Engineering Publications Limited. Ibrahmi, A., Bentaher, H., Hamza, E., Maalej, A., & Mouazen, A. M. (2015). Study the effect of tool geometry and operational conditions on mouldboard plough forces and energy requirement: Part 2. Experimental validation with soil bin test. Computers and Electronics in Agriculture, 117, 268e275. Li, S., Lobb, D. A., & Lindstrom, M. J. (2007). Tillage translocation and tillage erosion in cereal-based production in Manitoba, Canada. Soil and Tillage Research, 94(1), 164e182. Mak, J., Chen, Y., & Sadek, M. A. (2012). Determining parameters of a discrete element model for soiletool interaction. Soil and Tillage Research, 118, 117e122. Mari, I. A., Ji, C., Chandio, F. A., Arslan, C., Sattar, A., & Ahmad, F. (2015). Spatial distribution of soil forces on moldboard plough and draft requirement operated in silty-clay paddy field soil. Journal of Terramechanics, 60, 1e9. Raji, A. O. (1999). Discrete element modeliing of the deformation of bulk agricultural particles. PhD. University of Newcastle upon Tyne. Saunders, C. (2002). Optimising the performance of shallow, highspeed mouldboard ploughs. PhD. Cranfield University. Saunders, C., Godwin, R. J., & O’Dogherty, M. J. (2000). Prediction of soil forces acting on mouldboard ploughs. In Fourth international conference on soil dynamics. Adelaide, Australia. Ucgul, M., Fielke, J. M., & Saunders, C. (2014a). Three-dimensional discrete element modelling of tillage: Determination of a suitable contact model and parameters for a cohesionless soil. Biosystems Engineering, 121, 105e117. Ucgul, M., Fielke, J. M., & Saunders, C. (2014b). 3D DEM tillage simulation: Validation of a hysteretic spring (plastic) contact model for a sweep tool operating in a cohesionless soil. Soil and Tillage Research, 144, 220e227. Ucgul, M., Fielke, J. M., & Saunders, C. (2015). Three-dimensional discrete element modelling (DEM) of tillage: Accounting for soil cohesion and adhesion. Biosystems Engineering, 129, 298e306. Ucgul, M., Saunders, C., & Fielke, J. M. (2017a). Discrete element modelling of top soil burial using a full-scale mouldboard plough under field conditions. Biosystems Engineering, 160, 140e153.

70

b i o s y s t e m s e n g i n e e r i n g 1 9 0 ( 2 0 2 0 ) 5 8 e7 0

Ucgul, M., Saunders, C., & Fielke, J. M. (2017b). Discrete element modelling of tillage forces and soil movement of a one-third scale mouldboard plough. Biosystems Engineering, 155, 44e54. Vozka, P. (2007). Comparison of alternative tillage systems. School of Applied Science, Cranfield University (Masters by Research Thesis). Walton, O. (2006). Elastic-plastic contact model. Company Report. DEM Solutions. Walton, O. R., & Braun, R. L. (1986). Stress calculations for assemblies of inelastic spheres in uniform shear. Acta Mechanica, 63, 73e86.

Wang, X., Zhang, S., Pan, H., Zheng, Z., Huang, Y., & Zhu, R. (2019). Effect of soil particle size on soil-subsoiler interactions using the discrete element method simulations. Biosystems Engineering, 182, 138e150. Yan, Z., Wilkinson, S. K., Stitt, E. H., & Marigo, M. (2015). Discrete element modelling (DEM) input parameters: Understanding their impact on model predictions using statistical analysis. Computational Particle Mechanics, 2(3), 283e299. Zhang, L., Cai, Z., Wang, L., Zhang, R., & Liu, H. (2018). Coupled Eulerian-Lagrangian finite element method for simulating soil-tool interaction. Biosystems Engineering, 175, 96e105.