Prediction of Single Disc Seeding System Forces, Using a Semi-Analytical and Discrete Element Method (DEM)

Abstract

There is a rising interest amongst Australian farmers to use disc seeders due to their ability to operate in high residue conditions and at higher speeds, commonly in the range of 12 to 15 km h⁻¹. This paper reports on developing an analytical and discrete element method (DEM) force prediction model suited to a rotating flat disc blade operating at different sweep and tilt angles. To validate the models, field experiments were carried out with a flat disc blade at two tilt angles of 0 and 20° and four sweep angles of 6, 26, 45 and 90° in sandy soil. An analytical approach was developed following an experimental investigation that showed that only the forward portion of the disc blade is actively involved in generating soil failure, while the magnitude of this active portion of the soil-disc interface varied with sweep angle. The predicted active proportions correlated well with the experimental observations. As applying different sweep angles affects the direction of soil movement relative to the disc face, the directions of the friction and resultant forces at different sweep and tilt angles were determined. The equation of soil acceleration force was adapted to account for different sweep angles. Results showed that the predicted force fits relatively well with the measured data at 90, 45 and 26° sweep angle, while the low correlation between predicted and measured force at 6° sweep angle was due to the scrubbing reaction force not accounted for in the model. Results also showed that a better coefficient of determination (R² = 0.93) was obtained between DEM vs. test results compared to the analytical model predictions (R² = 0.86), particularly for predicting side forces. It was found from the study that both the developed analytical approach and DEM model enabled the prediction of soil forces at different sweep and tilt angles acting on a flat disc blade, which can assist in optimising disc design to lower the specific resistance.

1. Introduction

No-tillage (zero tillage) is widely adopted in Australia, where farms are large, and labour is expensive. Zero-tillage disc seeders can be operated at higher operation speeds and in heavier stubble residue conditions [1]. Accordingly, there has been increasing interest amongst Australian farmers in using disc seeders. However, the optimum disc angle settings of disc seeders and their impacts on soil forces are not well documented. This is particularly important to reduce the draught force requirements and hence the fuel consumption [2]. As a complement to field experiments which are time and resource-intensive and with seasonal limitations, suitable modelling methods can help speed up evaluation and improve the understanding of soil–disc interactions, to assist optimisation processes.

Analytical and numerical models are the two most common methods used to model soil–tool interactions. The interaction between soil and a concave disc was previously studied using analytical methods [3,4,5,6,7,8] based on wide blade theories [9,10,11]. Although these models provide some helpful information, they are limited as they assume that all of the disc–soil interfaces are involved in soil failure. This study developed a novel semi-analytical model considering that only the active part of the disc contributes to soil failure to predict the tillage forces under typical disc settings. Quantifying tillage forces accurately under different operation and geometry settings will contribute to improved disc seeder performance by determining optimum settings, which create less force without sacrificing the quality of the disc operation. The analytical model will also help designers and manufacturers who cannot afford to use expensive industry software to predict the tillage forces acting upon the disc seeder [11]. The basics of soil–disc blade interactions were experimentally reviewed under controlled soil conditions to help guide an improved basis of soil/disc blade modelling.

Many recent studies [12,13,14,15,16,17,18,19,20] have shown that the discrete element method (DEM) modelling can provide accurate soil force predictions when its parameters are accurately calibrated and applied. DEM was developed in the field of rock mechanics [21], and the technique is based on interacting particles that are able to simulate realistic soil bulk behaviour. Interactions between all particles are governed by contact models replicating physical laws and adjusted via a calibration process to suit specific soil conditions. Both tillage tool forces and soil movement can be predicted using DEM. The linear cohesion integrated hysteretic spring contact model can accurately simulate soil–tool interactions for tined sweep tillage tools [22], whereby contact forces between the particles are calculated using a hysteretic spring model, and cohesion is defined by adding a cohesive force between the particles. Although DEM is a valuable method, no comparison has been conducted to evaluate its ability to predict the tillage forces against an analytical model in the case of soil–disc interaction. Therefore, this study will provide valuable knowledge in comparing two different methods.

Logic of Approach

Initial observations and measurements of the soil loosening features created by a vertical, flat disc blade set at different sweep angles (The angle measured on a horizontal plane between the active face of the disc blade and the direction of travel) were conducted by the authors in an outdoor laboratory facility of the University of South Australia. The results (See Section 2.1) clearly showed that only the leading part of the disc blade was involved in generating the full extent of soil failure, while the trailing part contributed only to extra soil movement. This leading part of the disc blade could thus be described as the ‘active’ part of the blade, defined by its length at the soil surface.

From this initial work, the authors developed the following staged approach to model soil-to-flat disc blade interactions:

A.

The wide blade passive failure theory was first used to model a fixed circular blade operated at 90° to the direction of travel in a case study most closely related to wide blade theory, considering a wide circular blade with varying depth of cut across its width. The blade tilt angle (The angle measured between the active face of the disc blade and the vertical plane) was then the complementary angle to the blade rake angle. This basic model demonstrated the suitability of the wide blade theory for the case of a fixed circular blade operating in the range of 70 to 90° rake angle and validated a method of accounting for the variable depth of cut across a circular disc blade shape, as well as for the speed and bulldozing effects. The draught and vertical forces acting on the vertical blade were predicted within 6–19% and 1.5–14% of the measured data, respectively. This work was reported in [11,23].

B.

The above-fixed blade prediction model was extended to account for the additional effect of sweep angle with a focus on the 0 to 8° range combined with the effect of tilt angle in the 0 to 20° range. This combined range is directly relevant to the configurations of zero-tillage disc seeders using single disc openers. This model development stage introduced the concept of an active part of the blade involved in generating furrow soil failure and accounted for the soil surcharge accumulation developing at higher sweep angle values. These details are being reported in this paper.

C.

The impact of disc blade free-rotation was added to the stage 2 model with an analysis of the direction of soil flow impacting the direction of soil-tool frictional reactions.

D.

Finally, the effect of having a single-sided bevel on the rear face of the cutting edge was included to account for additional soil scrubbing reaction when operating at smaller sweep angles, which is often witnessed in field operations.

This paper describes the Stage 2 modelling in detail and validates the analytical force predictions with data obtained in a remoulded sandy-loam soil. It also compares the findings with DEM simulations of the same soil conditions, following the method reported by [22]. Comparing the two methods will help researchers and machinery manufacturers choose the suitable modelling approach.

2. Materials and Methods

2.1. Force Measurement Experiment

Experiments were performed in a continuous outdoor soil bin (tillage test track) located at the University of South Australia. During testing, the tillage test track had sandy loam soil (with no residue) at 7% moisture content (dry basis). The porosity of the soil was 39.5%. The particle size distribution of the soil is given in

Table 1. Particle size distribution of the test soil.

Sieve Size (mm)	Percentage Retained (%)
2.36	0.5
1.18	2.6
0.6	15.1
0.425	11.7
0.3	14.0
0.15	36.7
0.075	10.8
<0.075	8.6

.

Prior to testing, water was uniformly added and the soil was reconditioned via a loosening, levelling and multi-pass rolling process. A flat disc blade of 460 mm diameter, 6 mm thickness and wedge angle of 16° was used in the tests. All experiments were performed at a 67 mm disc operating depth within the range of typical seeding depth settings and at 2 km h⁻¹ forward speed, the lowest speed possible with the equipment, with the aim to approximate quasi-static conditions.

An experimental testing rig (Figure 1a) was used, which allowed the adjustment of the sweep angle (0 to 45° range) and tilt angle (−20 to +40° range). The draught, vertical, and lateral forces were simultaneously measured using a 3D force measuring frame with calibrated S-type load cells of 5 kN capacity (Figure 1b). The force measuring frame consisted of sets of 3 roller guided trays isolating each orthogonal force component, and its cross-sensitivities were calibrated under known static forces applied near the expected resultant force point of application on the disc blade. The forward speed was measured using a free-rolling side wheel (490 mm in diameter) equipped with a rotary encoder of 256 pulses/rev. The measured data were logged at a frequency of 100 Hz using LabView TM V 7.1 software.

Experiments were performed using four sweep angles (6, 26, 45 and 90°) and two tilt angles (0 and 20°) in a randomised complete block design with three (temporal) replications. Measurements were made over a distance of 40 m within the straight sections of the continuous, stadium-shaped track.

2.2. Active Disc Section Evaluation

To support the development of the analytical model, a series of tests were performed to determine the active section of a circular disc blade, which is involved in generating the full extent of soil failure. These exploratory experiments were performed in similar remoulded soil conditions and forward speed of 2 km h⁻¹. The experimental factors included three operating depths (40, 60 and 80 mm), two sweep angles (5 and 10°) and two tilt angles (0 and 20°), while contrasting the performance of a whole (Figure 2a) vs. a partial (Figure 2b) disc blade. The partial disc blade only had the forward 50% section of the disc blade engaged in the soil, while the rear section was cut out. The results of the active side evaluation are provided in the Supplementary Materials.

2.3. Analytical Model Development

2.3.1. Predicting the Active Chord Length, La

To quantify the active portion of the soil–disc interface, the theoretical approach outlined in [11,23] was used. The furrow width data measured by laser profile meter (e.g., Figure 3 and Figure 4) was used to determine the values of the variables La and m (the ratio of forward rupture distance f to operating depth d, both measured in situ) according to the geometrical relationship between the furrow boundary trace and the disc blade, at different sweep angles.

The chord length L was defined as the horizontal line intersecting the disc blade face at the undisturbed soil surface. In the modelling, the active chord length L_a was expected to extend beyond L/2 into the rear section of the disc blade. On the front section, the integration process was conducted over the variable Li (0 < Li < L/2) referenced from the disc blade centre (Li = 0 at full depth, d_max) and extending to a maximum L/2 when coinciding with zero operating depth. The integration process over the rear section was similarly conducted with the variable Li referenced from the disc blade centre to a maximum value (=L_a − L/2).

A geometrical analysis (Figure 3) defining L, L_i, L_a shows the geometrical relationship between furrow size (from disc blade point of entry B to full furrow width W_f. In the rear section of the disc (from L/2 toward L_a), soil forward rupture starts to reduce due to reducing of the depth of the disc. At the particular point E along the disc chord length located at a distance L_a from the disc blade point of entry B, the forward rupture distance f and the furrow width W_f reach their maximum value:

$$W_f=L_a\sin \beta +f \cos \beta$$

(1)

where $L_a$ is the active proportion of the disc chord (L) at the soil surface, and which is measured from the disc blade point of entry into the soil.

The reducing forward rupture distance f, further along the chord length, indicates that section of the disc blade does not contribute to soil failure.

The $L_a$ relationship with sweep angle [11] depicts a significant increase of $L_a$ over the range, from 0.5–0.55 L at the low sweep angle values encountered in zero-tillage disc seeders, and reaching 1.0 L at larger sweep angles (>70°). The value of $L_a$ can be used to predict the force at different sweep angles. For a flat disc opener set at a low sweep angle (e.g., 6°) as used on zero-till single disc seeders, the predicted relationship shows that the forward 55% portion of the soil–disc interface actively generates soil failure.

The integration process of passive soil failure force P is initiated from the disc centre at (d max) and thus is only valid for 0 < L_i < L/2. As shown in Figure 3, L_a is expected to lie in the range of L/2 to L; therefore, force equation is further modified conceptually as the sum of the segment two integrations, namely:

$$P={\displaystyle {\int }_0^{\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{P}_i{d}_l+{\displaystyle {\int }_0^{\left(L_a-\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.\right)}}{P}_i{d}_l$$

(2)

Letting k = L_a/L, Equation (2) can be re-written, as follows:

$$P={\displaystyle {\int }_0^{\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{P}_i{d}_l+{\displaystyle {\int }_0^{\left(k-\frac{1}{2}\right)L}}{P}_i{d}_l$$

(3)

For the passive soil failure reaction, P applied to the disc when the sweep angle β ≤ 90° [23] is outlined below:

$$\begin{array}{c}P=A_1\left[{\displaystyle {\int }_0^{\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{\left(\frac{b}{a}\sqrt{a^2-L_i{}^2}-\left(b-d\right)\right)}^{2}dl+{\displaystyle {\int }_0^{\left(k-\frac{1}{2}\right)L}}{\left(\frac{b}{a}\sqrt{a^2-L_i{}^2}-\left(b-d\right)\right)}^{2}dl\right](A_2+A_3+A_4\\ +A_5)\left[{\displaystyle {\int }_0^{\raisebox{1ex}{$L$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}\left(\frac{b}{a}\sqrt{a^2-L_i{}^2}-\left(b-d\right)\right)dl+{\displaystyle {\int }_0^{\left(k-\frac{1}{2}\right)L}}\left(\frac{b}{a}\sqrt{a^2-L_i{}^2}-\left(b-d\right)\right)dl\right]\end{array}$$

(4)

where:

$$A_1={\gamma }_i·N_{\gamma };$$

(5)

$$A_2=c·N_c;$$

(6)

$$A_3=c_a·N_{ca};$$

(7)

$$A_4=q·N_q;$$

(8)

$$A_5=\frac{{\gamma }_i}{g}·N_a{v}^2{\sin }^2\beta .$$

(9)

2.3.2. Validation of the Active Chord Length, La

To validate the L_a concept, an experiment was conducted in Tillage Test Track conditions at a range of operating depths, sweep and tilt angles, comparing a fixed whole disc blade with its assumed partial equivalent with the rear 50% of soil/tool interface removed (Figure 2). The three-dimensional soil reaction forces were measured in the Tillage Test Track environment using a replicated and randomised experiment. The results of measurements of draught, vertical and side force data for the whole and partial disc blades are shown in

Table 2. The measured forces of whole and partial disc blades at different angle settings.

Sweep Angle (°)	Tilt Angle (°)	Depth (mm)	Draught (N) Whole Partial		Vertical (N) Whole Partial		Side (N) Whole Partial
5	0	40	210	173	78	57	75	60
5	0	80	477	367	120	88	292	237
5	0	60	332	251	114	78	174	131
10	0	60	325	241	75	71	256	234
5	20	60	313	245	118	122	−19	−22
10	20	60	231	177	26	31	130	109

.

Table 2 shows that the measured forces were typically higher for the whole blade, namely 18 to 26% higher draught on average, 3 to 32% larger upward vertical force, and 4 to 9% greater side force. As the operating depth became shallower, draught and side force differences tended to become less significant. The results indicate that, while the forward half of the soil–disc interface generates most of the soil failure work, in practice, a slightly greater proportion is involved in creating the total soil/tool forces. Some differences in soil forces can be attributed to secondary soil/tool interaction involving the loosened soil being displaced sideways.

2.3.3. Accounting for the Tangential Frictional Reaction due to Sweep Angle

The wide blade passive failure theory as reported by [11,23] was used to model a fixed (non-rotating) circular blade set at 90° to the direction of travel and at two tilt angles of 0 and 20°. When the sweep angle β = 90°, the disc blade is normal to the direction of travel and thus the direction of soil–tool friction applies in the vertical plane only, combining with the effect of rake angle, as modelled in the wide blade soil passive failure theory, and influences the direction of the passive soil resultant force P. With a sweep angle β < 90°, the asymmetry of soil movement in the horizontal plane creates an additional tangential reaction that influences the direction of the resultant horizontal force component H_p (Figure 5). The forward movement of an angled disc blade was divided into orthogonal vector components that are perpendicular and tangential to the disc plane. The tangential component of blade movement generates a frictional reaction (${\mathsf{\delta }}_X$) at the soil/tool interface in the horizontal plane that is tangential to the disc face.

This component would be complementary to the traditional frictional reaction considered in the vertical plane (${\mathsf{\delta }}_Z,$ when β = 90°), which is mobilised by the perpendicular vector component of movement. The modified soil force modelling approach for a fixed blade at a sweep angle β < 90°, therefore, considers two elements of friction, namely δ_z in the vertical plane (Figure 5a) and δ_x in the horizontal plane (Figure 5b).

Observations with a fixed circular disc blade in the Tillage Test Track environment showed that the soil movement relative to the blade shifts from (i) a situation of relative movement occurring solely in the vertical plane at β = 90° to (ii) a situation where relative soil movement in the vertical plane reduces significantly with decreasing sweep angle, and (iii) a situation where movement in the horizontal plane predominates solely as the sweep angle approaches very low values. When β = $90°$, the following assumptions were therefore adopted δ_z = δ (friction fully mobilised in the vertical plane), and δ_x = 0° (friction not mobilised in the horizontal plane). When β is near 0°, the soil moves relative to the disc blade only in the horizontal direction and δ_x = δ (friction fully mobilised in the horizontal plane) and δ_z = 0° (friction not mobilised in the vertical plane). Similarly to [3], the transition in the δ_x and δ_z values from 90° sweep angle down to 0° was assumed to follow a reducing sine function for δ_z and an increasing cosine function for δ_x, as follows:

$${\mathsf{\delta }}_x=\delta \cos \beta$$

(10)

$${\mathsf{\delta }}_z=\delta \sin \beta$$

(11)

2.3.4. Accounting for the Effect of Velocity Differential

To account for the effect of velocity, the additional component ($\frac{{\gamma }_i}{g}·N_{a}d{v}^2$), as developed by [24], was added to the force equation. The inertial effect of forward velocity (v) on the disc forces assumed when β = 90° was reported in [23]. Additionally, it is necessary to account for the velocity differential occurring in the direction of the soil failure plane when inclined with a sweep angle. The following relationship was used for the passive soil failure velocity (v_p) of the disc:

$${\mathrm{v}}_P=v\sin \beta$$

(12)

At working speeds such as 20 km h⁻¹, these inertial effects have to be taken into account above a certain critical speed. This was suggested by [8] to be equal to $\sqrt{5gw}$, where g is the gravitational acceleration and w is the width of the tool.

In this study, the total force for the fixed blade was calculated over the sweep and tilt angle combinations using Equation (4) for the soil input parameters, including soil cohesion, adhesion, bulk density and soil surcharge. The effect of sweep angle on soil frictional reaction force was included in the partitioning of the soil passive soil failure force into its three force components, as illustrated in Figure 5 and calculated as follows:

$$H_p=P\sin (\alpha +{\delta }_z)$$

(13)

$$V_p=P\cos \left(\alpha +{\delta }_z\right)$$

(14)

$$S_p=H_p\cos (\beta +\mathsf{\delta }\mathrm{x})$$

(15)

$$D_p={\mathrm{H}}_{\mathrm{p}}\sin (\mathsf{\beta }+\mathsf{\delta }\mathrm{x})$$

(16)

The total draught force, $D_T$, in a bulldozing situation ($\mathsf{\beta }>45°)$ was calculated by adding a proportion (sin β) of the additional drag force $F$ and the horizontal component of the adhesion force R_a to the component D_p, as follows:

$$D_T=D_p +(F+R_a \cos \alpha )\sin \mathsf{\beta }$$

(17)

The total vertical force, $V_T$, in a bulldozing situation was estimated by combining the soil failure vertical force Vp, the soil body self-weight W₃, the interface force reaction H and the vertical component of the adhesion force R_a similarly to the methods expressed in [23].

The side force S_T was calculated as follows:

$$S_T=S_p +(F+R_a\cos \alpha )\cos \mathsf{\beta }$$

(18)

2.4. DEM Simulations

The DEM simulations were carried out using a DELL Precision T5810 Intel ^® Xeon CPU E5-2687W v4 @ 3.00 GHz computer with the software EDEM 2020^TM. The linear cohesion integrated hysteretic spring contact model allowed the particles to behave in a linear elastic manner up to predefined stress, and when the total stress on the contact area exceeded the predefined stress (=yield strength) in the model, the particles behaved as though undergoing plastic deformation. The cohesion between the particles was defined by adding a cohesion force to the normal contact forces. In the hysteretic spring contact model, the total normal (F_n) and tangential (F_t) forces were computed as

F_n = F_n^s + F_n^d

(19)

F_t = F_t^s + F_t^d

(20)

where F_n^s, F_n^d, F_t^s and F_t^d are the normal contact force, normal damping force, tangential contact force and tangential damping force, respectively. F_n^s was determined as per [25]:

$${\mathrm{F}}_{\mathrm{n}}^{\mathrm{s}}=-\{\begin{array}{llllll}{\mathrm{K}}_1·{\mathrm{U}}_{\mathrm{abn}}& & & & & \mathrm{loading}\\ {\mathrm{K}}_2·\left({\mathrm{U}}_{\mathrm{abn}}-{\mathrm{U}}_0\right)& & & & & \mathrm{unloading}/\mathrm{reloading}\\ 0& & & & & \mathrm{unloading}\end{array}$$

(21)

where U_abn is the normal component of the relative displacement, U_o is the residual overlap. K₁ and K₂ are the loading and unloading stiffnesses, respectively. As per [26,27], K₁ was computed as

K₁ = 5 r_eq min(Y_a, Y_b)

(22)

where Y is the yield strength, and r_eq is the equivalent radius and defined as where req is defined as [24]

1/r_eq= 1/r*_a+1/r*_b

(23)

where r* is the radius for the individual particles a and b. Following [26,27], K₂ was computed as

K₂ = K₁/e²

(24)

where e is the coefficient of restitution. The residual overlap was updated in each time step as

$${\mathrm{U}}_0=\{\begin{array}{llllll}{ \mathrm{U}}_{\mathrm{abn}}·\left(1-\frac{{\mathrm{K}}_1}{{\mathrm{K}}_2}\right)& & & & & \mathrm{loading}\\ {\mathrm{U}}_0& & & & & \mathrm{unloading}/\mathrm{reloading}\\ {\mathrm{U}}_{\mathrm{abn}}& & & & & \mathrm{unloading}\end{array}$$

(25)

The tangential contact force, F_t^s, was calculated as per [25] as

F_t^s = −n_k K₁ U_abt

(26)

where U_abt is the tangential component of the relative displacement. n_k is the stiffness factor defined as the ratio of tangential stiffness to normal loading stiffness. The normal and the tangential damping forces (F_n^d and F_t^d) were calculated using

F_n^d = −n_c (((4 m_eq K₁)/(1+(π/ln e)²)) Ů_abn)^−1/2

(27)

F_t^d = −(((4 m_eq n_k K₁)/(1+(π/ln e)²)) Ů_abt)^−1/2

(28)

where Ů_abn and Ů_abt are the normal and tangential components of the relative velocity, respectively. n_c is the damping factor that controls the amount of velocity-dependent damping. m_eq is the equivalent mass and is defined in [25] as

1/m_eq = 1/m*_a+ 1/m*_b

(29)

where m^* is the mass for the individual particles a and b. The total tangential force (F_t) was limited to the lesser of either the calculated tangential force or the sliding friction force. F_t was determined as

F_t = −min (n_k K₁ U_abt + F_t^d, µ F_n^s)

(30)

The magnitude of the moments caused by total tangential force (M) and the rolling resistance (M_r) were calculated following [28]:

M = r_con F_t

(31)

M_r = −µ_r F_n^s r_con λ_θ

(32)

where r_con is the perpendicular distance of the contact point from the centre of mass, μ_r is the coefficient of rolling friction, and λ_θ, is the unit vector of angular velocity at the contact point. The new position of a particle was computed by integrating Equations (33) and (34):

Ü = (F_n + F_t)/m*

(33)

R = (M + M_r)/I

(34)

The cohesion force was calculated as [25]:

F_c = ξ A_c

(35)

where ξ 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 A_c is the contact area. Hence, Equation (19) becomes

F_n = F_n^s + F_n^d + F_c

(36)

DEM is based on modelling the contacts between particles where many particles are used to represent a soil bulk. Hence, increasing the number of contacts increases the calculation time. In modelling soil, larger than actual spherical particle sizes are used to reduce the computation time (i.e., coarse graining). The DEM parameters were calibrated so that the larger than actual size DEM particles act as the actual soil. The DEM parameters used in the study are shown in

Table 3. EDEM parameters used for simulation.

Property	Value	Source
Density of sand particles (kg m−3)	2600	[29]
Density of steel (kg m−3)	7861	[30]
Shear modulus of soil (Pa)	5 × 107	[31]
Shear modulus of steel (Pa)	7.9 × 1010	[30]
Poisson’s ratio of soil	0.3	[32]
Poisson’s ratio of steel	0.3	[33]
Coefficient of restitution of soil-soil	0.1	Calibrated
Coefficient of friction of soil-soil	0.7	[34]
Coefficient of friction of soil-steel	0.5	[34]
Coefficient of rolling friction of soil-soil	0.28	Calibrated
Coefficient of rolling friction of soil-steel	0.05	[34]
Cohesive energy density between soil-soil (N m−2)	9000	Measured using DST *
Adhesive energy density between soil-tool (N m−2)	3700	Measured using DST *

* DST = Direct shear test.

, comprising a combination of measured and literature-based parameters. The DEM parameters of the coefficient of rolling friction and coefficient of restitution between soil and soil (a value less than 0.3 for a compressible soil) were calibrated by matching simulation results to the actual measurement result of an angle of repose test.

In order to measure the angle of repose, a soil sample was placed in a pipe (100 mm diameter and 300 mm long). After that, the pipe was lifted upward (at around 500 ms⁻¹). Soil flowed into a cylindrical tray (200 mm diameter and 22.5 mm long) until the soil overflowed and formed a pile. When at rest, the soil’s angle of repose flowed on the surface and was measured using a digital level. The test procedure was also simulated in DEM.

By using the calibrated parameters (given in Table 3), an angle of repose of 29° (measured in the test) was achieved in the simulation using a trial and error process (with an error margin of ±1° i.e., −28.2° measured in DEM)

In order to run the simulations, a virtual soil bin of 4000 mm long × 1500 mm wide × 300 mm deep was generated. A nominal 40 mm particle size was used with the particles generated based on the measured particle size distribution of the sandy loam soil. A total of 359,188 particles were used in the simulations.

The DEM bulk density was set to match the bulk density used in the experiments. To do so, particles were compressed using an upper physical plane until the desired bulk density was achieved.

The disc models were created using SolidWorks^TM software and then imported into the EDEM environment. After the depth and the speed of the disc models were set, simulations were carried out (Figure 6). The tillage forces were averaged over the stable section of the force vs. displacement graphs.

3. Results

The relationships between force components and sweep angles of a fixed disc are presented in Figure 7 and Figure 8 for the blade at two tilt angles and operating at a low speed (2 km h⁻¹). The analytically predicted draught and vertical forces correlate well with the measured data over the different sweep angles (Figure 9a). Draught force expectedly increased with a greater sweep angle due to increasing the active part of the disc and associated furrow disturbance (Figure 7a and Figure 8a). The rate of increase was lower at the 20° tilt angle (Figure 8a), but, overall, was well below that shown in the literature for concave discs (e.g., [3]). This can be explained by the lower amount of soil displacement induced by a flat disc blade compared to concave discs of harrows and ploughs used in previous literature, which inverts a large amount of soil. In line with experimental observations, the modelled effect of bulldozing (=accumulation of soil body ahead of the blade) becomes increasingly significant as the sweep angle increases from 45 to 90°.

The upward vertical force increases with a greater sweep angle, explained by the proportionality of the vertical reaction to increasing soil passive failure reaction force with an increasing sweep angle (Figure 7b and Figure 8b). This strongly contrasts with the decreasing relationships observed by [3,35], who found a decreasing upward vertical reaction by increasing the sweep angle with concave discs. Their findings can be attributed to the effect of disc concavity, inducing (i) a much stronger push upward at low sweep angles from the compacting scrubbing reaction on the convex face of the blade and (ii) the beneficial effect of lower apparent rake angle on the concave face as the sweep angle increases. The modelled effect of bulldozing on the vertical force also becomes increasingly significant as the sweep angle increases from 45 to 90° while being slightly lower than that for the draught force. It was also found that the analytical model yields downward vertical forces (instead of the upward vertical forces measured) for a 20° tilt angle (Figure 8b).

The results for the side force show an expected rise to a maximum at approximately 30 to 40° sweep angle and a decrease beyond that (Figure 7c and Figure 8c). Contrary to the literature data for concave discs, the side force does not reverse in direction at low sweep angles due to the limited scrubbing reaction occurring at the bevel edge angle, but overall, was well below that shown in the literature; e.g., [3] with concave discs. This can be explained by the lower extent of soil displacement induced by a flat disc blade compared to concave discs of harrows and ploughs used in previous literature, which invert a large amount of soil.

In line with experimental observations, the modelled effect of bulldozing becomes increasingly significant as the sweep angle increases from 45 to 90° (Figure 7c and Figure 8c). The side force becomes nearly equal to the draught force when $\beta +{\mathsf{\delta }}_X=45°$. The analytical model significantly overpredicts the side force at the lower sweep angles, as the effect of scrubbing on the rear of the blade is not taken into account at this stage of the model. Overall, the draught and vertical forces acting on the vertical blade at 2 km h⁻¹ speed were predicted for the four experimental sweep angles within 17% to 27% of the measured values. As detailed above, the side force was significantly over-predicted at sweep angles below 45° due to scrubbing. The model accurately reflected the beneficial effect of a 20° tilt on reducing the draught force and on the side force at the larger sweep angles (Figure 7 vs. Figure 8). The vertical force reaction was significantly under-predicted at lower sweep angles as the scrubbing reaction at the bevel edge was not accounted for in the model. This scrubbing would increase the vertical reaction force. Prediction for the side force fell within 6.7% of the measured values at sweep angles of 45° and above. There was a significant overprediction of side force at lower sweep angles.

The DEM predicted draught, vertical and side forces of the vertical disc showed a good correlation with the measured data over the range of different sweep angles even though the analytical solution performs better for some angles. Compared to the analytical results, a better coefficient of determination was obtained using DEM simulations (R² = 0.93 vs. R² = 0.86) (Figure 9). A better side force prediction was obtained using DEM simulations, while vertical forces were slightly over-predicted for 6 and 25° sweep angles. The vertical force over-prediction at sweep angles below 45° can be attributed to the larger than actual size particles used in the simulations. In the draught and side force predictions, the DEM particles contact the disc’s large surface, while only a few are in contact with the disc’s surface in the vertical force predictions. An extensive scaling effect study, which was not considered in this study, might also help to determine contact parameters more precisely when the particle sizes are scaled up. Additionally, using clump particles (a combination of multiple spheres to obtain a more realistic particle shape) might also increase the accuracy of the simulations. In addition, in this study, the calibration procedure was based on an angle of repose test; performing a second test for the calibration procedure, such as a penetration test, might help to improve the accuracy of the DEM parameters and hence predicted vertical force prediction.

4. Conclusions

In this study, a fixed circular blade was tested in remoulded sandy soil at a range of sweep and tilt angles. An analytical soil reaction force model previously reported [11,23] was modified to suit the case of a flat circular disc blade (fixed in rotation) operating at various sweep and tilt angles suited to disc seeder configurations. It was found that the analytical model developed can predict draught forces very well, whereas it does not accurately predict the vertical and side forces for sweep angles of 6°. To improve the results of the analytical model future work will need to account for the scrubbing reaction at the rear bevel edge and to account for disc blade-free rotation.

Such model can then be used in sensitivity analyses to identify optimum points of operation including specific draught resistance, tilt and sweep angles and bevel edge geometry under contrasting soil conditions. Opportunities for further work are also identified, including modelling the association of disc blade and depth gauge wheel, often used on many zero-till disc seeders.

DEM simulation of the soil and fixed flat disc blade interaction was also carried out. The DEM results over-predicted the draught force at high sweep angles for a 0° tilt angle and overpredicted the draught force for both low and high sweep angles at a tilt angle of 20°. The DEM gave reasonable predictions of upward vertical force. However, the DEM predictions gave much higher side forces than that measured at sweep angles of 6 and 25° for both of the tilt angles of 0 and 20°. This was despite the DEM model accounting for the scrubbing reaction that occurs on the rear disc face at low sweep angles. Further work is needed to improve the model and adapt it to a freely rotating flat disc blade.