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

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

<sup>2</sup> Faculty of Science and Engineering, Southern Cross University, Lismore, NSW 2084, Australia

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

**Abstract:** Disc seeders are commonly used in no-till farming systems, and their performance evaluation generally rely on expensive and time-consuming field experiments. Mathematical models can help speed up force-related evaluations and improve the understanding of soil-disc interactions, to assist the performance optimisation processes. Previous analytical force prediction models of disc blades have not accounted for the free rotation aspect of the disc blade. This paper develops an analytical force prediction model from the wide blade failure theory adapted to suit rotating flat disc blades operating at different sweep and tilt angles and compares predictions with Discrete Element Method (DEM) simulations. To validate the two models, experiments were performed on a remoulded sandy soil condition using a rotating flat disc set at two tilt angles of 0◦ and 20◦, and four sweep angles of 6, 26, 45 and 90◦ the 3-dimensional force components of draught, vertical and side forces were measured. Results showed a higher coefficient of determination (R<sup>2</sup> = 0.95) was obtained with analytical model predictions compared to DEM predictions (R2 = 0.85) for their agreement with the test results. It was found that both the developed analytical approach and the DEM model can be used to predict tillage forces at different sweep and tilt angles acting on a rotating flat disc blade.

**Keywords:** flat disc; analytical force prediction model; discrete element method (DEM); soil-tool interaction

#### **1. Introduction**

During the last two decades, Australian no-till farmers have increasingly adopted disc seeders due to their characteristics of low soil disturbance and low soil throw, which facilitates working at higher speeds and operating at narrow rows spacings with high amounts of residue [1,2]. However, the disc seeder soil-reaction forces, particularly the draught force, must be optimised to minimise power requirements and fuel consumption [1] at higher speeds. Designing and evaluating tillage tools are generally based on time and resource-intensive field tests, which can only be undertaken at certain times of the year. If the interaction between soil and disc can be modelled, a large majority of the field tests can be avoided.

Currently, there is no published analytical force prediction model for a flat rotating disc blade with a low sweep angle (3–8◦) and tilt angle (0–20◦), as used in Australian no-till disc seeders. The analytical force prediction models developed by [3–7] modelled spherical/conical disc blades with large sweep angles (15–80◦) and tilt angles (up to 30◦), as normally fitted on implements used for soil inversion and residue incorporation. Limited work has been conducted using DEM for modelling soil-disc interaction by [8,9].

In past work by the authors [10], the interaction between soil and disc was modelled using an analytical method and a discrete element method (DEM) model. For the analytical method, tests were first conducted in an outdoor soil bin which revealed that only the leading part of the disc (active part) was involved in generating soil failure. Based on this

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

Academic Editor: Tao Cui

Received: 15 December 2022 Revised: 6 January 2023 Accepted: 10 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/).

<sup>1</sup> STEM, University of South Australia, Adelaide, SA 5000, Australia

information and using the wide blade passive failure theory [11,12], a fixed/non-rotating circular blade, operating at 90◦ to the direction of travel at two tilt angles of 0 and 20◦, was modelled to prove the suitability of the wide blade passive failure theory for modelling a flat circular disc [13,14]. In a second step, the authors extended the fixed blade model to account for the effect of different sweep (0 to 8◦) and tilt angles (0 to 20◦) relevant to the configurations of zero-tillage disc seeders as used on single disc openers. Although the model provided some helpful information at that stage, it was limited as it assumed that the disc did not rotate.

As an alternative approach, the interaction between soil and a fixed disc blade was also modelled using DEM, which was proven as an effective method when its parameters were accurately calibrated [15–19]. In this study, the modelling work of [10] is extended to predict the draught, vertical and side forces reactions applied onto a flat rotating disc blade. Both analytical and DEM predictions are compared to the experimental results.

#### **2. Materials and Methods**

#### *2.1. Experimental Work*

The tillage forces (draught, vertical and side forces) acting on a circular disc were measured in Tillage Test Track (TTT) at the University of South Australia. The TTT is a continuous outdoor soil bin. It has two straights of 50 m length joining together by two curves by 50 m diameter. The test soil was placed between rails (2.5 m wide by 0.3 m deep). In the TTT, a tractor tows two trolleys, each capable of tillage testing [20]. A hydraulic cylinder on the test frame was used to adjust the frame height to obtain the desired operating depth (Figure 1). The TTT had sandy loam soil (85% sand, 3% silt and 12% clay with no residue), and during testing, it had a 7% moisture content (dry basis). The soil was reconstituted by tilling the full width, levelling and consolidating it using a 2 t troller.

**Figure 1.** Tillage test track (TTT) rig adapted from [20].

To perform the tests, a flat disc blade of 460 mm diameter, 6 mm thickness and wedge angle of 16◦ was used in the tests. Experiments were performed at 67 mm operating depth (a typical seeding depth in Australia) and 2 km h−<sup>1</sup> forward speed. The operating speed of the disc seeders widely used in Australia is usually in the range of 12 to 15 km h−1. However, this was undertaken to establish an analytical model based on the wide blade theory developed for quasi-static conditions. The lowest practical speed for testing of 2 km h−<sup>1</sup> was therefore selected to approximate quasi-static conditions and matched in the development of the analytical model. A test frame was developed that allowed setting different sweep (0 to 45◦) and tilt (−20 to 40◦) angles (Figure 2). The tillage forces were measured using a 3D dynamometer frame with 5 kN capacity S-type load cells. The data was logged at a frequency of 100 Hz using LabView TM V 7.1 software. Tests were performed for four sweep angles (6, 26, 45 and 90◦) and two tilt angles (0 and 20◦) using a randomised complete block design with three replications. Measurements were conducted over 40 m within the straight sections of the TTT (where a steady state condition was reached).

**Figure 2.** (**a**) Force measurement frame fitted with adjustable disc mounting frame and (**b**) the definition of sweep and rake angles.

#### *2.2. Analytical Model Development*

The rotation of the disc was observed during testing to affect the direction of soil flow. Thus, the analytical model accounted for the change in soil flow by changing the direction of action of friction and resultant force. This affected the magnitude and direction of its components of draught, vertical and side forces.

#### 2.2.1. To Determine the Active Part of the Disc

The active portion of the disc involved in soil failure was calculated to calculate the soil reaction forces onto the disc, as explained in [13]. Due to the difference in the *m* ratio (the ratio of forward rupture distance f to operating depth d, both measured in situ) between a fixed and free rotating disc at different sweep angles, the value of the active part of the disc, *La* needed to be determined.

The active part of the disc (*La*) and ( *La <sup>L</sup>* ) ratio was determined and is shown in Figure 3. As shown in Figure 3, the calculated ( *La <sup>L</sup>* ) ratios at different sweep angles were similar to the measured values.

**Figure 3.** Comparison of calculated and measured ( *La <sup>L</sup>* ) ratio at α = 90◦ (α is the rake angle).

2.2.2. Calculation of the Force

Once the active part of the disc, *La*, was determined, the method explained in [13,14] was used to calculate the forces at different sweep and tilt angles. The calculated forces at rake angles of α = 90◦ and 70◦ at 2 km h−<sup>1</sup> forward speed were compared with the measured forces with the difference shown in Table 1.

**Table 1.** Analytical force predictions compared to measured results.


The large difference between the predicted and measured force at a sweep angle of 6◦ can be explained by the rear side scrubbing reaction force, which was not accounted for by the model used in this work. Thus, the wide blade theory predicts well the force of a rotating disc blade only at the larger sweep angles.

Since the rotation of the disc affects the direction of soil flow and the direction of action of friction, the friction components of δ*<sup>x</sup>* and δ*<sup>z</sup>* were calculated and used to determine the force components of the rotating disc's draught, vertical and side forces.

#### 2.2.3. Effect of Disc Rotation on the Direction of Soil Reaction Forces

With a free-rotating disc, the rotation of the disc changes the soil flow direction and thus affects the direction of friction at different sweep angles. The direction of friction alters the direction of the resultant force, P. To determine the effect of rotation on the direction of force, P, the direction of soil flow must be determined. The direction of overall soil movement over the disc face depends on the result of pressures applied on the disc surface due to sweep angle, disc rotation and the resistance of soil to compressive failure. It was assumed soil moves over the disc face as a rigid block. The direction of friction can be assumed as opposed to the direction of soil movement. Thus, the elemental friction forces act in the opposite direction to the velocity vectors [21]. To determine the overall soil movement direction, the relative movement of the soil to the disc, which is a function of both the forward and rotational speed of the disc, was determined. In order to determine the average soil-disc relative velocity, the soil–disc relative velocity at any point on the disc surface was calculated. A kinematic analysis of the disc [21] showed a considerable variation in the relative velocity over the soil-disc contact area. Therefore, the point of action of the resultant force on the soil-disc contact area was assumed to be the point at which the relative velocity is equal to the average relative velocity. That means the friction vector at this point is representative of the average friction vectors that apply to the soil-disc contact area.

2.2.4. Method to Determine the Point of Action of force P on the Soil-Disc Contact Zone

In order to determine the point of action of the resultant force P on the soil-disc interface plane, the point of action of force P on the Z (vertical plane) and X (horizontal plane) axis was determined. The joining of these two points is the point of action of force P on the disc-soil contact area.

Method to Determine the Point of Action of Force P on a Vertical Plane

To determine the point of action of force P on a vertical plane, the value and point of action of any individual force acting on this plane were determined. For this, the percentages of the required force for soil failure in each element (e.g., the ratio of required force to overcome soil cohesion, surcharge, etc., to the total force, P) at different sweep angles were required. All applied forces on the active side of the disc–soil contact area are shown in Figure 4. It was assumed that the cohesion force, Pc, adhesion force, Pca and surcharge, Pq, act at half of the operating depth and the force related to soil weight and speed act at 2/3 of the depth [11].

**Figure 4.** Horizontal and vertical applied force on the disc-soil contact area of the active side of the disc to determine the point of action.

The percentages of required force for each element were calculated from the force calculation equation [13]. Therefore, by knowing the value and point of action of any individual force, the point of action of a resultant elementary force (Pi) on a vertical plane was determined, as shown in Equation (1).

$$\begin{array}{c} \frac{(P\_c + P\_{ci} + P\_q)P\_i}{2} \ast d\_i + \frac{2(P\_\gamma + P\_a)P\_i}{3} \ast d\_i = P\_i \ast b\_{zi} \\ b\_{zi} = \frac{3(P\_i + P\_{ci} + P\_q) + 4(P\_\gamma + P\_a)}{6} d\_i \end{array} \tag{1}$$

Since depth (di) is variable over the disc chord, the point of action of the resultant force P on a vertical plane (*bz*), was calculated using Equations (2) and (3) and the integration equation described in [6] by replacing the r, d and L in the equation.

$$b\_z = \frac{3(P\_c + P\_{c1} + P\_{q)} + 4(P\_\gamma + P\_a)}{6} \ast \frac{P\_1 d\_1 + P\_2 d\_2 + \dots \cdot \dots P\_n d\_n}{P} \tag{2}$$

$$b\_{z} = \frac{\Im(P\_{\varepsilon} + P\_{ax} + P\_{q)+} 4\left(P\_{\gamma} + P\_{\bar{x}\bar{\beta}}A\_{1} \int\_{0}^{L\_{\bar{x}}} \left[\left(\sqrt{a^{2} - L\_{x}^{2}} - (b-d)\right)^{2} \cdot d\_{i}\right] dl + \left(A\_{2+}A\_{3} + A\_{4+}A\_{5}\right) \int\_{0}^{L\_{\bar{x}}} \left[\left(\sqrt{a^{2} - L\_{x}^{2}} - (b-d)\right) \cdot d\_{i}\right] dl}{6P} \tag{3}$$

where *<sup>A</sup>*<sup>1</sup> <sup>=</sup> *<sup>γ</sup><sup>i</sup>* ·*Nγ*, . *<sup>A</sup>*<sup>2</sup> <sup>=</sup> *<sup>c</sup>*·*Nc*, *<sup>A</sup>*<sup>3</sup> <sup>=</sup> *ca*·*Nca* , *<sup>A</sup>*<sup>4</sup> <sup>=</sup> *<sup>q</sup>*·*Nq*, *<sup>A</sup>*<sup>5</sup> <sup>=</sup> *<sup>γ</sup><sup>i</sup> <sup>g</sup>* ·*Na* sin *<sup>β</sup>v*2.

As *La* > L/2, the above calculation was performed first by substituting *La* with L/2 and then with (*La* − L/2) as per Equation (4).

$$\begin{aligned} b\_{2} = \frac{3\left(P\_{i} + P\_{a} + P\_{q} + \right)4\left(P\_{\gamma} + P\_{a}\right)\left[A\_{1}\int\_{0}^{\frac{L}{2}}\left[\left(\sqrt{a^{2} - L\_{a}^{2}} - (b - d)\right)^{2} \cdot d\_{i} + \int\_{0}^{\frac{L}{2} - L\_{a}}\left(\sqrt{a^{2} - L\_{a}^{2}} - (b - d)\right)^{2} \cdot d\_{i}\right] dt + \int\_{0}^{\frac{L}{2}}\left(\sqrt{a^{2} - L\_{a}^{2}} - (b - d)\right)^{2} \cdot d\_{i}\right] dt}{\left(A\_{2} + A\_{3} + A\_{4} \int\_{0}^{\frac{L}{2}}\right)\left[\left(\sqrt{a^{2} - L\_{a}^{2}} - (b - d)\right) \cdot d\_{i} + \int\_{0}^{\frac{L}{2} - L\_{a}}\left(\sqrt{a^{2} - L\_{a}^{2}} - (b - d)\right) d\_{i}\right]\right)} \end{aligned} \tag{4}$$

To Determine the Point of Action of Force P on a Horizontal Plane

It was assumed that *P*1, *P*2, ... .*Pn* are the forces acting at a distance *ax*1, *ax*2, ... .*axn* from the Z axis. The Z axis was chosen as a reference axis. The point of action of the different applied forces (*P*1, *P*2, ... .*Pn* ) on a horizontal plane (on X-axis direction, distance *ax*) were determined by the summation of force times distance from the reference axis (Z) and dividing it by the total force P applied on the assumed plane, as shown in Equations (5)–(8).

$$P.a\_x = \sum\_{0}^{L\_x} (P\_i a\_{xi})\tag{5}$$

$$a\_{\!\!\!0} = \frac{P\_1 a\_1 + P\_2 a\_2 + \cdots + P\_n a\_n}{P} \tag{6}$$

$$\sigma\_x = \frac{A\_1 \int\_0^{L\_a} \left[ \left( \sqrt{a^2 - L\_a^2} - (b - d) \right)^2 . L\_i \right] dl + \left( A\_{2+} A\_3 + A\_4' + A\_5 \right) \int\_0^{L\_a} \left[ \left( \sqrt{a^2 - L\_a^2} - (b - d) \right) . L\_i \right] dl}{P} \tag{7}$$

As *La* > L/2, the above calculation was performed first by substituting *La* by L/2 and then by (*La* − L/2) as per Equation (8).

$$a\_{\lambda} = \frac{A\_1 \int\_0^{\frac{1}{2}} \left[ \left( \sqrt{a^2 - L\_a^2} - (b - d) \right)^2 \cdot d + \int\_0^{\frac{1}{2}} \left( \sqrt{a^2 - L\_a^2} - (b - d) \right)^2 \cdot d\_i \right] dl +}{P}$$

$$\left[ (A\_{2 +} A\_3 + A\_{4 + A\_5} \int\_0^{\frac{1}{2}} \left[ \left( \sqrt{a^2 - L\_a^2} - (b - d) \right) \cdot d\_i + \int\_0^{\frac{1}{2}} \left( \sqrt{a^2 - L\_a^2} - (b - d) \right) d\_i \right] dl \right]$$

Thus, point B, the joining point of action of forces P in X and Z directions (Figure 4), is the point of action of force P.

#### 2.2.5. Method to Determine the Direction of Friction

Once the point of action of the resultant force P was determined (it was assumed an average resultant velocity act at this point), to determine the direction of soil-disc friction, it was required to calculate the resultant velocity at point B, as it was assumed that the elemental friction forces act in the opposite direction to the velocity vectors [21]. Therefore, a moving system of coordinates OXZ was used to study the kinematics of the disc [21]. By rotating the axis OX, through a sweep angle, β on a horizontal plane and rotating the axis OZ, through a tilt angle, α<sup>t</sup> on a vertical plane, the new coordinate system of OX1Z1 will be formed (Figure 5). The axis OX is along with the direction of travel, and axis OX1 and OZ1 are in the plane of the disc cutting edge. The component of the resultant speed (Vr) in the X direction at point B is:

$$V\_{\rm rx} = V - V\_l \cos(\theta 0^\circ - \theta) \tag{9}$$

$$V\_t = V\lambda \left(r - (d - b\_z)\right) / rV\_t\\ \text{remains constant, but } V \text{ at OX1 becomes } V \cos \beta \tag{10}$$

$$V\_{rx1} = V\cos\theta - (V\_{l}\sin\theta) \tag{11}$$

$$V\_{rz1} = V\_l \cos \theta \tag{12}$$

$$\theta = \tan^{-1} \frac{AO}{AB} \tag{13}$$

$$V\_r = \sqrt{V\_{rx1}^2 + V\_{rz1}^2} \tag{14}$$

where *V* = forward speed at the direction of travel, Vt = tangential speed at point B, *Vr* = resultant speed (absolute velocity), *λ* = speed ratio from the experimental data.

By substituting the value for *V*, β and θ, the resultant velocity was calculated. The measured data for a rotating disc at different sweep angles for λ (speed ratio) was used to calculate the tangential speed (Equation (10)).

By calculating the direction of average soil movement relative to the disc at different sweep angles, the angles *λ<sup>x</sup>* and *λ<sup>z</sup>* were determined (Figure 6) as follows:

$$\cos \lambda\_x = \frac{V\_{r\ge 1}}{v\_r} \tag{15}$$

$$
\cos \lambda\_z = \frac{V\_{rz1}}{v\_r} \tag{16}
$$

**Figure 6.** The direction of resultant soil movement friction.

The direction of force, P, that was assumed perpendicular to the disc face will be changed by the action of friction. In the wide blade theory, it is assumed that soil moves only upward on the disc face; therefore, friction only changes the direction of force P in a vertical plane, but for a rotating flat disc with a sweep angle, the soil moves with an angle relative to the horizontal and vertical planes (Figures 5 and 6). Friction has components in both vertical and horizontal directions (Figure 7). Therefore, it changes the direction of force, P, on both vertical (Z-axis) and horizontal planes (X-axis). The effect of friction on the orientation of force, P relative to a vertical plane, was determined by calculating the angle that force P makes with a vertical plane perpendicular to the direction of travel (δz) (Figure 7a). The effect of friction on the orientation of force P relative to the horizontal plane was determined by calculating the angle that force P makes with a vertical plane parallel to the direction of travel (δx) (Figure 7b). When the disc sweep angle is 90◦ the situation is the same as a wide blade and soil moves vertically without a horizontal movement, so δ<sup>z</sup> = δ and δ<sup>x</sup> = λ*<sup>z</sup>* =0 (Figure 6). When the disc sweep angle is zero, soil moves only on a horizontal plane, and no vertical movement occurs, so δ<sup>x</sup> = δ, λ*<sup>z</sup>* = 90◦ and λ*<sup>x</sup>* = 0. Hence, the sine function was used to show the λ*<sup>z</sup>* and λ*<sup>x</sup>* variation as follows:

$$
\delta\_{\mathbf{x}} = \delta \ast \sin \lambda\_{\mathbf{z}} \text{ Or } \delta\_{\mathbf{x}} = \delta \ast \cos \lambda\_{\mathbf{x}} \tag{17}
$$

$$
\delta\_z = \delta \ast \cos \lambda\_z \tag{18}
$$

By calculating the friction angle components on a vertical and horizontal plane, the direction of force P and their components (draught, vertical and side forces) at different sweep angles were determined.

#### 2.2.6. Method to Determine the 3D Force Component of the Free-Rolling Disc

Figure 7 shows how the force P was partitioned into its draught, vertical and side force components at different sweep angles.

#### *2.3. DEM Simulations*

In this study, 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. A linear cohesion integrated hysteretic spring contact model was suggested by [16] to model the soil–circular disc opener interactions. The hysteretic spring contact model accounts for the plastic deformation behaviour in the contact mechanics equations. Compressible materials such as soil can be modelled using the hysteretic spring contact model resulting in particles behaving in a linear elastic manner up to a predefined stress (Figure 8). When the total stress on the contact area exceeds the predefined compressive stress in the model, the particles behave as though they are undergoing plastic deformation [22]. On the other hand, the linear cohesion model allows users to add a cohesive force to the normal force direction with the ability to model soil cohesion.

**Figure 8.** Hysteretic spring force-displacement relationship [14].

The hysteretic spring contact model calculates the total normal (*Fn*) and tangential (*Ft*) forces as;

$$F\_{\text{II}} = F\_{\text{II}}{}^s + F\_{\text{II}}{}^d \tag{19}$$

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

*Fn <sup>s</sup>* was found as per [22];

$$\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}$$

According to [23,24], *K1 and K2* were calculated, as;

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

$$K\_2 = K\_1 \text{\textquotedblleft} \text{\textquotedblright} \tag{23}$$

where *req* is [22],

$$1/r\_{eq} = 1/r\_{\
u}^{\*} + 1/r\_{\
u}^{\*} \tag{24}$$

The residual overlap was updated in each time step as;

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

*Ft s* , *Fn d,* and *Ft <sup>d</sup>* were calculated as per [22];

$$F\_t^{\;s} = -n\_k \, K\_1 \, \mathcal{U}\_{\text{alt}} \tag{26}$$

$$F\_{\rm II}^{\rm d} = -n\_{\rm c} \left( ((4 \text{ } m\_{\rm cq} \text{ K}\_1)/(1 + (\pi/\ln \, e)^2)) \, \text{ $\vec{\mathcal{U}}\_{\rm abm}$ } \right)^{-1/2} \tag{27}$$

$$F\_{\rm I}^{\rm d} = -(((4\ \,\mathrm{m\_{tq}}\ \,\mathrm{n\_{k}}\ \,\mathrm{K\_{1}})/(1+(\pi/\ln\,e)^{2}))\,\mathrm{U}\_{\rm abt})^{-1/2}\tag{28}$$

where *meq* is the equivalent mass and was defined in [22] as;

$$1\,\mathrm{/m}\_{cq} = 1\,\mathrm{/m\stackrel{\circ}{\_{a}}} + 1\,\mathrm{/m\stackrel{\circ}{\_{b}}} \tag{29}$$

*Ft* was calculated using the following equation;

$$F\_l = -\min\left(n\_k \,\, \mathcal{K}\_I \,\, \mathcal{U}\_{abt} + F\_l^{\;\,d}, \,\mu \, F\_n^{\;\,s}\right) \tag{30}$$

*M* and *Mr* were computed as suggested by [25]:

$$\mathcal{M} = r\_{con} \, F\_t \tag{31}$$

$$M\_r = -\mu\_r \, F\_n{}^s \, r\_{\rm conv} \, \lambda\_\theta \tag{32}$$

The new position of a particle was updated by integrating the following equations;

$$
\dot{\Omega} = (F\_m + F\_t) / m^\* \tag{33}
$$

$$R = (M + M\_r)/I \tag{34}$$

The cohesion force was calculated as [22];

$$F\_{\mathfrak{c}} = \mathfrak{J}^{\mathfrak{c}} \mathcal{A}\_{\mathfrak{c}} \tag{35}$$

As a result, Equation (18) becomes;

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

Using actual particulate sizes and shapes in 3D DEM simulations increases computation time (due to the high number of contacts) and computation costs. Therefore, particle sizes used in 3D DEM are generally selected larger than the actual particle sizes. In modelling soil, larger than actual spherical particle sizes are used to reduce the computation time. Hence, appropriate mechanical properties must be defined via calibration for the bulk behaviour of these larger particles. The DEM parameters used in this study are shown in Table 2. The DEM parameters of the coefficient of rolling friction and coefficient of restitution between soil and soil (a value less than 0.3 for compressible soil) were calibrated by matching simulation results to the actual measurement result of an angle of repose test.

**Table 2.** DEM parameters used for simulation.


\* DST = direct shear test.

In order to measure the angle of repose, a soil sample was placed in a pipe (100 mm diameter and 300 mm long). The pipe was then lifted upward (at around 500 mm s<sup>−</sup>1). Soil flowed into a cylindrical tray (200 mm diameter and 22.5 mm high) until the soil overflowed and formed a pile. When at rest, the soil's angle of repose was measured using a digital level. The same test procedure was simulated in DEM. By using the calibrated parameters (Table 2), a static angle of repose of 29◦ (measured) was achieved in the simulation using a trial and error process (with an error margin of ±1◦ with 28.2◦ measured in DEM). The results of the trial and error process are shown in Table 3.

**Table 3.** Calibration results angle of repose simulation.


\* selected values.

DEM simulations were undertaken using a virtual soil bin with 4000 mm length, 1500 mm width and 300 mm depth. A nominal 40 mm particle size (spherical) with the same particle size distribution of the test soil was used (359,188 particles) (Figure 9). The DEM bulk density was set to match the bulk density used in the experiments of 1860 kg m−<sup>3</sup> (dry basis). This was achieved by compressing the DEM particles using an upper physical plane until the bulk density of the virtual soil media reached the bulk density used in the field tests. The disc models were created using SolidWorksTM software and imported as a STEP file into the EDEM software. An operating depth of 67 mm and speed of 2 km h−<sup>1</sup> were used. As the circular disc is a passive-driven tool (as the soil's force drives its rotation on the disc), the rotational velocities measured during the TTT experiments (shown in Table 4) were used in the DEM simulations. In the future, the DEM simulation could be coupled with multibody dynamic (MBD) software to find the disc's rotational speed. From the simulations, the tillage forces were measured over the steady state section of the force vs. displacement graphs.

**Figure 9.** DEM simulation of the soil-disc interaction.

**Table 4.** Rotational speed measured in the TTT and used in the DEM simulations.


#### **3. Results**

The relationships between the measured and predicted (both analytically and DEM) draught, vertical and side forces and sweep angles of the rotating disc are presented in Figures 10 and 11 for 0 and 20◦ tilt angles.

**Figure 10.** Comparison of measured analytically predicted and DEM predicted forces for a 0◦ tilt angle. The measured values include error bars for 1 standard deviation.

The test results showed that the accumulated soil in front of the disc at a higher sweep angle adds an additional force to the draught force while the side force first increased with the increase of sweep angle (up until 26◦) and then reduced steeply as the disc pushes the soil forward rather than to the side. Also, the cosine function is reduced by increasing the angle. The vertical force decreased until the 26◦ sweep angle and then significantly increased with the increase of the sweep angle. Tilt angle has a major effect on soil reaction force, especially on the vertical force, as the tilted disc soil reaction force has downward components that help the disc penetrate the soil. Due to the close relationship between the rake angle and tilt angle in the disc, the tilt angle affects force partitioning in 3 dimensions. Overall, by tilting the disc 20◦, the applied draught force was reduced by 19% at 6◦ and

40% at 90◦ sweep angle while the vertical force was decreased by 25% at 6◦ and 82% at a 90◦ sweep angle. Increasing the tilt angle to 20◦ also reduced the side forces by 70% at 6◦ and 38% at 90◦ sweep angles.

The analytically predicted draught, vertical and side force results showed a similar trend to the test results over the different sweep angles, but a large difference between the draught force component at β = 90◦ and the side force component at β = 6◦ was observed. The large difference between predicted and measured side force data at β = 6◦ can be explained by scrubbing reactions which were not considered in this work. In comparison, the difference between the measured and predicted draught force β = 90◦ can be attributed to the high amount of soil accumulated in front of the disc, which could not be accurately accounted for in the analytical model.

The DEM predicted draught, vertical and side forces of the vertical disc generally showed a good correlation with the measured data over the range of different sweep angles; however, DEM could not predict the draught forces at a 20◦ tilt angle well. When compared to the analytical results, a better coefficient of determination was obtained using analytical results simulations (R2 = 0.95) than the DEM results (R2 = 0.85) (Figure 12). The lower coefficient of determination for the DEM simulation results may be able to be improved by further refining the DEM calibrated parameters to better suit the larger than actual particles and using DEM-multibody dynamic coupling.

**Figure 12.** The whole of the data set correlations between the measured reaction forces and (**a**) analytical predictions and (**b**) DEM predictions (n = 24).

#### **4. Conclusions**

In this study, an analytical soil-force model initially developed for a fixed flat disc blade in [10] was extended to predict the forces for a rotating flat disc blade over a range of sweep and tilt angles. Its predictions were compared with experimental measurements, and results showed that this extended analytical model could accurately predict the soil reaction forces (R<sup>2</sup> = 0.95). The rotating flat disc blade interactions with a similar soil condition were carried out in simulated DEM space by assigning the disc blade rotational speed as a constant value. At this stage, the DEM could not predict the draught forces as accurately as the analytical model at the higher tilt angles. However, the benefits of DEM simulation include soil movement predictions and therefore can deliver more holistic results for optimising disc blade performance. Using multibody coupling to allow reactive simulations rather than assigning a fixed rotational velocity may improve accuracy and would also allow the DEM simulations to be extended beyond the quasi-static work conducted to date. Therefore, future research should focus on the multibody dynamic-DEM coupling to improve and extend the analysis. Beyond the current work conducted in low-moisture sandy-loam soil, it would be beneficial to undertake a further study into sticky soil conditions, a common challenge for most disc seeders.

**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., C.S. 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.
