*2.2. The Discrete Element Model of Rotary Blade-Soil–Straw*

#### 2.2.1. Soil and Straw Model Parameters

The discrete element parameters consist of the intrinsic parameters (material density, Poisson's ratio, shear modulus), material contact parameters (coefficient of restitution, coefficient of static friction, dynamic friction coefficient), and contact model parameters. The parameters employed in this article are taken from the reference [21]. The relevant parameters are shown in Table 2.

The increase in the number of particles in the software will lead to a decline in the geometric level of the simulation calculation speed [22]. In addition, too small particles will lead to an uncertain impact on the model under the micro force. Therefore, the size of the soil model is generally larger than the real soil particles in the simulation calculation. Without affecting the simulation test, the soil particle radius selected in this paper is 5 mm, 10 mm, 15 mm, and 30 mm, as shown in Figure 3.

**Figure 3.** Soil particles with different sizes.

In this paper, the rotary tillage blade mainly covers the straws scattered on the ground during operation and will not crush and cut the straws. For the purpose of improving the simulation efficiency, the straws established in this paper are also simplified to a certain extent. It is mainly composed of 10 particles with a radius of 5 mm. The center distance between adjacent particles is 5 mm. The total length of the straw model is 55 mm. The simulation model of straw is shown in Figure 4.

**Figure 4.** Straw model.

**Table 2.** Material physical properties parameters.


#### 2.2.2. Contact Model Settings

Edem software 2018 integrates many options to apply to different material types of contact. We selected Hertz-Mindlin with Bonding along with the characteristics that the soil particles are difficult to reconnect after fracture in this paper. At the same time, the most commonly used contact theoretical model Hertz-Mindlin (no slip), is selected for ordinary contact calculation.

#### (a) Hertz-Mindlin (no slip)

In this model, the normal force component model is based on the Hertzian contact theory. The tangential force component model is based on the work of Mindlin-Deresiewicz. Both normal and tangential forces have damping components, where the damping coefficient is related to the coefficient of restitution. The tangential friction force follows the Coulomb friction law, and the rolling friction adopts the research results of Sakaguchi. The relevant calculation principles are described in the literature [23–28].

The normal force *Fn* is expressed as follows:

$$F\_{\rm II} = \frac{4}{3} E \ast \sqrt{R^\*} \delta\_n^{\frac{3}{2}} \tag{8}$$

$$\frac{1}{E^\*} = \frac{(I - v\_i^2)}{E\_i} + \frac{(1 - v\_j^2)}{E\_j} \tag{9}$$

$$\frac{1}{R^\*} = \frac{1}{R\_i} + \frac{1}{R\_j} \tag{10}$$

where *Fn* is the normal force, with the unit of N, *E*\* is equivalent to Young's modulus, with the unit of Pa, *R*\* is the equivalent radius, with the unit of m, *δ<sup>n</sup>* is the amount of normal overlap, with the unit of m, *vi,j* is Poisson's ratio, *Ei,j* is Particle Young's modulus, with the unit of Pa, *Ri,j* is Particle radius, with the unit of m.

The normal damping force *Fn <sup>d</sup>* is expressed as follows:

$$F\_n^d = -2\sqrt{\frac{5}{6}} \beta \sqrt{S\_d m^\*} v\_n^{\rightarrow} \tag{11}$$

$$\beta = \frac{-\ln e}{\sqrt{\ln^2 e + \pi\_n^2}}\tag{12}$$

$$S\_n = 2E^\* \sqrt{R^\* \delta\_n} \tag{13}$$

where *Fn <sup>d</sup>* is Normal damping force, with the unit of N, *β* is viscous damping coefficient, *e* is coefficient of restitution, *Sn \_*is normal stiffness, with the unit of N/m, *m*\* is equivalent

mass, with the unit of g, *mi,j* is Particle mass, with the unit of g, *v* → *rel <sup>n</sup>* is Normal relative velocity, with the unit of m/s.

The expression for the tangential force *Ft* is as follows:

$$F\_t = -S\_t \delta\_t \tag{14}$$

$$S\_t = 8G^\* \sqrt{R^\* \delta} \tag{15}$$

where *Ft* is the normal force, N,*St* is tangential stiffness, with the unit of N/m, *δ<sup>t</sup>* is tangential overlap, with the unit of m, *G*\* is equivalent shear modulus, with the unit of Pa.

The tangential damping force *Ft <sup>d</sup>* is expressed as follows:

$$F\_t^d = -2\sqrt{\frac{5}{6}} \beta \sqrt{S\_t m^\*} \upsilon\_t^{\stackrel{\rightarrow}{rel}} \tag{16}$$

where *Ft <sup>d</sup>* is angential damping force, with the unit of N, *St* is tangential stiffness, with the unit of N/m, *v* → *rel <sup>t</sup>* is Tangential relative velocity, with the unit of m/s.

The total tangential force is limited by the Coulomb friction force. When the tangential force reaches the Coulomb friction force, that is, *Ft* > *μsFn*, relative sliding occurs between the particles, where *<sup>s</sup>* is the static friction coefficient. To simplify the contact model, The model does not take into account the sliding friction of the contact model. The rolling friction is expressed in terms of the contact torque of the particle model, and the expression is as follows:

$$
\pi\_i = -\mu\_r F\_n R\_i \omega\_i \tag{17}
$$

where *τ<sup>i</sup>* is torque, with the unit of N·m, *μ<sup>r</sup>* is coefficient of rolling friction, *Ri* is the distance from the contact point to the centroid, with the unit of m, *ωi*—angular velocity of particles at the contact point, with the unit of r/min.

### (b) Hertz-Mindlin with Bonding

Hertz-Mindlin with Bonding contact model bonds particles together. Since this model includes Hertz-Mindlin (no slip), the default Hertz-Mindlin (no slip) contact model needs to be deleted from the software Creator Tree list when generating the bond.

Hertz-Mindlin with Bonding contact model uses finite-sized "bonding" bonds to bond particles. The bond can withstand the normal and tangential motion between the contacting particles to generate normal and tangential shear stress until the bond breaks at the critical maximum stress. The particles interact as rigid spheres, and no new bonds are generated, which is in line with the soil stress mechanical behavior characteristics of remaining loose after crushing [29–31].

In this model, the first generated particles interact with each other according to the Hertz-Mindlin contact model, and the connection starts after the set bond generation time point. After the particles are connected, the force and torque between the particles are both set to 0, according to the following Stepwise adjustment at each time step:

$$
\delta F\_n = -\upsilon\_n S\_n A \delta\_t \tag{18}
$$

$$
\delta F\_t = -\upsilon\_t S\_t A \delta\_t \tag{19}
$$

$$
\delta M\_{\rm ll} = -\omega\_{\rm n} S\_{\rm t} I \delta\_{\rm t} \tag{20}
$$

$$
\delta M\_t = -\omega\_t S\_n \frac{I}{2} \delta\_t \tag{21}
$$

$$A = \pi R\_B^2 \tag{22}$$

$$J = \frac{1}{2}\pi R\_B^4\tag{23}$$

where *Fn,t* is normal and tangential force, with the unit of N, *vn,t* is normal and tangential velocity, with the unit of m/s, *Sn,t* is normal and tangential stiffness, with the unit of N/m, *A* is the cross-sectional area of bonding bond, with the unit of m2, *δ<sup>t</sup>* is time step, with the unit of s, *Mn,t* is normal and tangential moments, with the unit of N·m, *ωn,t* is normal and tangential angular velocities, with the unit of r/min, *J* is the moment of inertia, with the unit of m4, *RB* is the bond-bond radius.

When the normal and tangential stresses exceed preset critical values, the bonded bonds will break, the particles will interact as rigid spheres, and no new bonds will be formed. Spawn bonds in this model can create cohesive forces when particles are generated, even if the particles are not in contact, so the contact radius should be set larger than the radius of the particle body.

Table 3 lists the parameters adopted in this study. In this article, we use Edem's dynamic generation method to generate soil particle beds layer by layer, and the model is shown in Figure 5.

The particle contact model can quantify the interaction between discrete elements within a given contact range. Adding a contact model to the particle model is a key step in the test. Appropriate contact model to accelerate simulation test.


**Table 3.** Contact parameters.

**Figure 5.** The discrete element model.

2.2.3. Initial Conditions for Simulation

In this simulation test, in order to ensure that the tillage depth was equal to 200 mm, the rotation center was set to the corresponding position into the soil according to the radius of rotation of each blade, as shown in Figure 6.

**Figure 6.** Initial position of the rotary blade.

According to the research objectives of this article, combined with the existing research experience and actual working conditions, the main indicators were determined as the power (P), the rate of crushed soil (RC), and the rate of straw burial (RB). Figure 1b shows that the value of B is fixed when R and β are certain, so the factor B is removed. The initial conditions of the experiment are shown in Table 4, and the factor level design is shown in Table 5.


**Table 4.** Initial conditions.

#### **Table 5.** Factor levels.


In this study, for convenience of calculation, the ratio of the number of bond-bond breaks per unit volume to the number of bond bonds at the beginning is used as the index of soil breakage. The sampling position is set 500 mm from the center of rotation after rotary tiller operation, and 1000 mm is set by grid-grid function after EDEM treatment. 1000 mm × 500 mm × 250 mm grid solver calculates the bonding bond between particles as a problem to be solved.

The value of power (P) is selected from the data within 2 s of the operation stabilization stage. Each blade is tested three times, and the average value is obtained, which is directly exported after EDEM software post-processing.

The crushed soil rate (RC) is measured by the ratio of bond breakage between particles in the tillage area. In this study, for the convenience of calculation, the ratio of the number of soil bond breakage per unit volume to the number of bonded bonds at the initial time was used as the indicator of crushed soil rate. The sampling position was set at 500 mm from the center of rotation after the rotary tillage knife operation, and the grid solver of 1000 mm × 500 mm × 250 mm was set using the grid function of EDEM post-processing, and the particle-to-particle bond was calculated as the problem to be solved.

Straw burial rate (RB) is the percentage of straw residual mass on the ground posttillage to the pre-tillage. In the post-processing of the EDEM software, L0 is the straw layer. The grid solver was divided into four layers (L1, L2, L3, L4) in the vertical direction, and the percentage of the total straw mass from L1 to L4 to the initial L0 mass was calculated as the straw burial rate.
