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

**Ali Khosravani 1, Jacky M. A. Desbiolles 1, John M. Fielke 1, Mustafa Ucgul 2,\* and Chris Saunders <sup>1</sup>**


**\*** Correspondence: mustafa.ucgul@scu.edu.au

**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<sup>−</sup>1. 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<sup>2</sup> = 0.93) was obtained between DEM vs. test results compared to the analytical model predictions (R<sup>2</sup> = 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.

**Keywords:** disc seeder; disc blade; discrete element method (DEM); force prediction; semi-analytical model; sandy soil

#### **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–8] based on wide blade theories [9–11]. Although these models provide some helpful information, they are limited as they assume that

**Citation:** Khosravani, A.; Desbiolles, J.M.A.; Fielke, J.M.; Ucgul, M.; Saunders, C. Prediction of Single Disc Seeding System Forces, Using a Semi-Analytical and Discrete Element Method (DEM). *Agriculture* **2023**, *13*, 206. https://doi.org/ 10.3390/agriculture13010206

Academic Editor: Massimo Cecchini

Received: 14 December 2022 Revised: 7 January 2023 Accepted: 11 January 2023 Published: 13 January 2023

**Copyright:** © 2023 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

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–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:


failure and accounted for the soil surcharge accumulation developing at higher sweep angle values. These details are being reported in this paper.


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.

**Table 1.** Particle size distribution of the test soil.


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−<sup>1</sup> 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.

**Figure 1.** Equipment used for disc blade experiments (**a**) disc blade adjustment rig and (**b**) threedimensional force measuring frame.

#### *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−1. 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.

**Figure 2.** Example soil failure from full (**a**) and partial (50%) (**b**) disc blade soil engagement, while fixed in rotation.

#### *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., Figures 3 and 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 *La* 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, *dmax*) 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 (=*La* − *L*/2).

A geometrical analysis (Figure 3) defining *L*, *Li*, *La* shows the geometrical relationship between furrow size (from disc blade point of entry B to full furrow width *Wf*. In the rear section of the disc (from *L*/2 toward *La*), 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 *La* from the disc blade point of entry B, the forward rupture distance f and the furrow width Wf reach their maximum value:

$$\mathcal{W}\_f = L\_a \sin \beta + f \cos \beta \tag{1}$$

where *La* 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 *La* relationship with sweep angle [11] depicts a significant increase of *La* 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 *La* 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 < *Li* < *L*/2. As shown in Figure 3, *La* 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 = \int\_0^{1/2} P\_l d\_l + \int\_0^{(L\_d - l/2)} P\_l d\_l \tag{2}$$

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

$$P = \int\_0^{l/2} P\_i d\_l + \int\_0^{(k-\frac{1}{2})L} P\_i d\_l \tag{3}$$

**Figure 3.** Plan view outlining the geometrical relationships from the record of furrow boundary trace relative to the disc blade at a sweep angle β and represented by its chord length L at the soil surface.

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

$$\begin{split} P = A\_1 \Big[ \int\_0^{l/2} \left( \frac{b}{a} \sqrt{a^2 - L\_l^2} - (b - d) \right)^2 dl + \int\_0^{(k - \frac{1}{2})L} \left( \frac{b}{a} \sqrt{a^2 - L\_l^2} - (b - d) \right)^2 dl \Big] (A\_2 + A\_3 + A\_4) \\ &+ A\_5) \left[ \int\_0^{l/2} \left( \frac{b}{a} \sqrt{a^2 - L\_l^2} - (b - d) \right) dl + \int\_0^{(k - \frac{1}{2})L} \left( \frac{b}{a} \sqrt{a^2 - L\_l^2} - (b - d) \right) dl \right] \\ &\qquad \times \int\_0^{l/2} \left( \frac{b}{a} \sqrt{a^2 - L\_l^2} - (b - d) \right) dl + \int\_0^{(k - \frac{1}{2})L} \left( \frac{b}{a} \sqrt{a^2 - L\_l^2} - (b - d) \right) dl \right] \end{split} (4)$$

where:

$$A\_1 = \gamma\_{\bar{\imath}} \cdot N\_{\bar{\imath}'} \tag{5}$$

$$A\_2 = c \cdot N\_c;\tag{6}$$

$$A\_3 = c\_a \cdot N\_{ca};\tag{7}$$

$$A\_4 = q \cdot \mathcal{N}\_{q'} \tag{8}$$

$$A\_5 = \frac{\gamma\_i}{\mathcal{S}} \cdot \mathcal{N}\_d \upsilon^2 \sin^2 \beta. \tag{9}$$

**Figure 4.** Example average soil furrow profile measured at 6◦ sweep angle of a full disc blade. Blue line indicates the furrow area and green line indicates the x-coordinate for drawing.

#### 2.3.2. Validation of the Active Chord Length, La

To validate the *La* 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.

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


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 *Hp* (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 (δ*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 (δ*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 δ<sup>z</sup> in the vertical plane (Figure 5a) and δ<sup>x</sup> 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 δ<sup>z</sup> = δ (friction fully mobilised in the vertical plane), and δ<sup>x</sup> = 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 δ<sup>x</sup> = δ (friction fully mobilised in the horizontal plane) and δ<sup>z</sup> = 0◦ (friction not mobilised in the vertical plane). Similarly to [3], the transition in the δ<sup>x</sup> and δ<sup>z</sup> values from 90◦ sweep angle down to 0◦ was assumed to follow a reducing sine function for δ<sup>z</sup> and an increasing cosine function for δx, as follows:

$$
\delta\_x = \delta \cos \beta \tag{10}
$$

$$
\delta\_z = \delta \sin \beta \tag{11}
$$

**Figure 5.** Partitioning of the passive soil reaction force P on a tilted disc blade at a sweep angle β; (**a**) side view (along the plane of the disc) and (**b**) top view.

#### 2.3.4. Accounting for the Effect of Velocity Differential

To account for the effect of velocity, the additional component ( *<sup>γ</sup><sup>i</sup> <sup>g</sup>* ·*Nadv*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 (*vp*) of the disc:

$$\mathbf{v}\_P = \upsilon \sin \beta \tag{12}$$

At working speeds such as 20 km h−1, these inertial effects have to be taken into account above a certain critical speed. This was suggested by [8] to be equal to 5*gw***,** 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\_{\mathcal{P}} = P \sin(\mathfrak{a} + \delta\_{\mathfrak{z}}) \tag{13}$$

$$V\_p = P \cos \left(\alpha + \delta\_z\right) \tag{14}$$

$$S\_p = H\_p \cos\left(\beta + \delta \mathbf{x}\right) \tag{15}$$

$$D\_p = \mathbb{H}\_{\mathbb{P}} \sin(\mathfrak{\beta} + \mathfrak{\delta} \mathfrak{x}) \tag{16}$$

The total draught force, *DT*, in a bulldozing situation (β > 45◦) was calculated by adding a proportion (sin *β*) of the additional drag force *F* and the horizontal component of the adhesion force *Ra* to the component *Dp*, as follows:

$$D\tau = D\_p \, \, + \, \left( F + R\_d \, \, \cos \, n \right) \sin \beta \, \, \tag{17}$$

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

The side force *ST* was calculated as follows:

$$S\_T = S\_{p^\*} + (F + R\_d \cos \alpha) \cos \beta \tag{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 2020TM. 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 (*Fn*) and tangential (*Ft*) forces were computed as

$$F\_n = F\_n{}^s + F\_n{}^d \tag{19}$$

$$F\_l = F\_l{}^s + F\_l{}^d \tag{20}$$

where *Fn s* , *Fn <sup>d</sup>*, *Ft <sup>s</sup>* and *Ft <sup>d</sup>* are the normal contact force, normal damping force, tangential contact force and tangential damping force, respectively. Fn <sup>s</sup> was determined as per [25]:

$$\mathbf{F\_n^s} = -\begin{cases} \mathbf{K\_1} \cdot \mathbf{U\_{abn}} & \text{loading} \\ \mathbf{K\_2} \cdot (\mathbf{U\_{abn}} - \mathbf{U\_0}) & \text{unloading/reloading} \\ 0 & \text{unloading} \end{cases} \tag{21}$$

where *Uabn* is the normal component of the relative displacement, *Uo* is the residual overlap. *K*<sup>1</sup> and *K*<sup>2</sup> are the loading and unloading stiffnesses, respectively. As per [26,27], *K*<sup>1</sup> was computed as

$$\mathbb{K}\_1 = \mathbf{5} \text{ } r\_{\text{eq}} \text{ } min(\text{Y}\_{a\text{\prime}}, \text{Y}\_b) \tag{22}$$

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

$$1/\text{r}\_{\text{eq}} = 1/\text{r}\_{\text{a}}^{\*} + 1/\text{r}\_{\text{b}}^{\*} \tag{23}$$

where *r\** is the radius for the individual particles *a* and *b*. Following [26,27], *K*<sup>2</sup> was computed as

$$K\_2 = K\_1 / \epsilon^2 \tag{24}$$

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

$$\mathbf{U}\_{0} = \begin{cases} \mathbf{U}\_{\text{abn}} \cdot \left( 1 - \frac{\mathbf{K}\_{1}}{\mathbf{K}\_{2}} \right) & \text{loading} \\ \mathbf{U}\_{0} & \text{unloading/reloading} \\ \mathbf{U}\_{\text{abn}} & \text{unloading} \end{cases} \tag{25}$$

The tangential contact force, *Ft s* , was calculated as per [25] as

$$F\_I{}^s = -n\_k \ K\_1 \downarrow I\_{\text{alt}} \tag{26}$$

where *Uabt* is the tangential component of the relative displacement. *nk* is the stiffness factor defined as the ratio of tangential stiffness to normal loading stiffness. The normal and the tangential damping forces (*Fn <sup>d</sup>* and *Ft <sup>d</sup>*) were calculated using

$$\mathbf{F\_n}^{\rm d} = -\mathbf{n\_c} \left( \left( \left( 4 \text{ m}\_{\rm eq} \,\mathrm{K}\_1 \right) / \left( 1 + \left( \pi / \ln \,\mathrm{e} \right)^2 \right) \right) \,\stackrel{\circ}{\mathrm{U}\_{\rm abn}} \right)^{-1/2} \tag{27}$$

$$\mathbf{F\_t^d} = -\left( \left( \left( 4 \text{ m}\_{\text{eq}} \text{ n}\_{\text{k}} \text{ K}\_1 \right) / \left( 1 + \left( \pi / \ln \text{ e} \right)^2 \right) \right) \mathbf{\hat{U}\_{\text{abt}}} \right)^{-1/2} \tag{28}$$

where *U˚ abn* and *U˚ abt* are the normal and tangential components of the relative velocity, respectively. *nc* is the damping factor that controls the amount of velocity-dependent damping. *meq* is the equivalent mass and is defined in [25] as

$$1/m\_{a\eta} = 1/m^\*\_{a} + 1/m^\*\_{b} \tag{29}$$

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

$$F\_t = -m \text{in } (n\_k \,\, \mathcal{K}\_1 \,\, \mathcal{U}\_{\text{alt}} + F\_t \, ^d \,, \, \mu \, F\_n \, ^s) \tag{30}$$

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

$$\mathcal{M} = r\_{\rm com} \, F\_{\rm t} \tag{31}$$

$$\mathcal{M}\_{\mathbb{F}} = -\mu\_{\mathbb{F}} \, F\_{\mathbb{H}} \, ^{s} \, r\_{\text{con}} \, \lambda\_{\theta} \tag{32}$$

where *rcon* is the perpendicular distance of the contact point from the centre of mass, *μ<sup>r</sup>* 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):

$$
\overrightarrow{CI} = (F\_{\text{fl}} + F\_{\text{l}})/m^\* \tag{33}
$$

$$R = (M + M\_{\rm r}) / I \tag{34}$$

The cohesion force was calculated as [25]:

$$F\_{\mathfrak{C}} = \mathfrak{J}^{\mathfrak{c}} \mathcal{A}\_{\mathfrak{C}} \tag{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 *Ac* is the contact area. Hence, Equation (19) becomes

$$F\_{\rm ll} = F\_{\rm n}{}^{s} + F\_{\rm n}{}^{d} + F\_{\rm c} \tag{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, 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−1). 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)

**Table 3.** EDEM parameters used for simulation.


\* DST = Direct shear test.

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 SolidWorksTM 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.

**Figure 6.** DEM simulation of the soil–disc interaction in a cohesionless medium.

#### **3. Results**

The relationships between force components and sweep angles of a fixed disc are presented in Figures 7 and 8 for the blade at two tilt angles and operating at a low speed (2 km h−1). 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 (Figures 7a and 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 (Figures 7b and 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 (Figures 7c and 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◦ (Figures 7c and 8c). The side force becomes nearly equal to the draught force when *β* + δ*<sup>X</sup>* = 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−<sup>1</sup> 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 (R2 = 0.93 vs. R2 = 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.

**Figure 7.** Comparison of measured, analytically predicted, and DEM predicted forces for 0◦ tilt angle (**a**) draught forces, (**b**) vertical forces, and (**c**) side forces. Error bars for measured results are plus and minus standard error.

**Figure 8.** Comparison of measured, analytically predicted and DEM predicted forces for 20◦ tilt angle. (**a**) draught forces, (**b**) vertical forces, and (**c**) side forces.

**Figure 9.** Whole data set of coefficients of determination between measured reaction forces and (**a**) analytical predictions or (**b**) DEM predictions (*n* = 24). The red dots represent the measured force vs. analytically predicted force (as (x,y)), and the blue dots represent the measured force vs. DEM predicted force (as (x,y)). On the other hand, the black line represents the perfect prediction (R<sup>2</sup> = 1) which was added to visually show how closely the models predicted the measured data.

#### **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.

**Supplementary Materials:** The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/agriculture13010206/s1, Table S1: Measured force and furrow area of whole/partial disc blades at different settings. Figure S1: Measured draught force of (Full) whole/partial disc blades at different settings. Figure S2: Measured vertical force of whole/partial disc at different settings. Figure S3: Measured side force of whole/partial (cut) disc at different settings. Figure S4: Comparison of furrow area of full and partial disc.

**Author Contributions:** Conceptualisation, A.K., J.M.A.D. and J.M.F.; methodology, A.K. and J.M.A.D.; Software (DEM) assistance: M.U. and A.K.; validation, A.K.; formal analysis, A.K. and J.M.A.D.; writing—original draft preparation, A.K.; writing—review and editing, M.U., C.S., J.M.A.D. and J.M.F.; supervision, J.M.A.D. and J.M.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Grains Research and Development Corporation (GRDC) (Project Number: USA00005) and the University of South Australia (UniSA).

**Institutional Review Board Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available upon request from the first author A.K.

**Acknowledgments:** The authors gratefully acknowledge the Grains Research and Development Corporation (GRDC) and the University of South Australia (UniSA) in co-funding for this PhD research.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Nomenclature**



#### **References**


**Disclaimer/Publisher's Note:** The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.
