Simulation of Mouldboard Plough Soil Cutting Based on Smooth Particle Hydrodynamics Method and FEM–SPH Coupling Method

Abstract

In the field of agricultural machinery, various empirical field tests are performed to measure the tillage force for precision tillage. However, the field test performance is costly and time-consuming, and there are many constraints on weather and field soil conditions; the utilization of simulation studies is required to overcome these shortcomings. As a result, the SPH method and the coupled FEM-SPH method are used in this paper to investigate the mouldboard plough–soil interaction. In this paper, the finite element software LS-DYNA was used to build the SPH model and the FEM-SPH coupling model of soil cutting, as well as to investigate the change in cutting resistance during the soil cutting process. The simulation results are compared with those of the experiments, and the curves of the simulation and experiment are in good agreement, which verifies the reliability of the model. The validated simulation model was used to investigate the effects of the cutting speed, depth of cut, inclination angle, and lifting angle of the mouldboard plough on cutting resistance. The simulation studies show that the SPH model takes 5 h and 2 min to compute, while the FEM-SPH coupled model takes 38 min; obviously, the computational efficiency of the FEM-SPH coupled model is higher. The relative errors between the SPH model and the experiment are 2.17% and 3.65%, respectively. The relative errors between the FEM-SPH coupled model and the experiment are 5.96% and 10.67%, respectively. Obviously, the SPH model has a higher computational accuracy. The average cutting resistances predicted by the SPH model and the FEM-SPH coupled model, respectively, were 349.48 N and 306.25 N; these resistances are useful for precision tillage. The cutting resistance increases with the increase in cutting speed and is quadratic; the cutting resistance increases with the increase in cutting depth and is quadratic; the horizontal cutting resistance and the combined cutting resistance increase with the increase in inclination angle, while the vertical cutting resistance remains essentially constant with the increase in inclination angle; the horizontal cutting resistance and combined cutting resistance increase as the lifting angle increases, while the vertical cutting resistance decreases as the lifting angle increases.

1. Introduction

The mouldboard plough is still today’s most widely used soil tillage tool [1,2]. Since Roman times, the mouldboard plough has been widely used as the dominant tillage tool for the following purposes: soil turning, burying litter and crop residues, laying the foundation for seedbeds, and making the soil loose and airy [3]. However, ploughing with a plough is a highly energy-consuming process [4,5], and the high energy consumption of agricultural machinery results in the emission of large amounts of exhaust gases, which cause pollution [6]. Precision agriculture is a technology-enabled farming management approach that observes, measures, and analyzes the needs of individual fields. Precision tillage involves measuring within-field soil strength variations and making applications to the field accordingly rather than applying an equal amount of cutting force during tillage operations in each field. This approach significantly reduces tillage energy consumption [7]. As a result, it is necessary to investigate the mouldboard soil tillage process and its mechanical properties to precisely till the soil and thus reduce energy consumption.

In recent years, many researchers have investigated the soil–tool interaction and the mechanical properties. Given that tillage force is primarily related to soil properties, tillage speed, tillage depth, and tool geometry [8], scholars have investigated the effects of soil properties [9,10], tillage speed [11,12], tillage depth [13,14], and tool geometry [15,16] on tillage force using field experiments and other methods. Some researchers have also developed tillage prediction models [17,18]. However, developing a fast and accurate tillage prediction model remains difficult because of differences in soil types, working conditions, and other factors.

Experimental and analytical methods are the two approaches to studying soil–tool interactions [19]. The experimental method is intuitive and reliable. However, the experimental method requires significant labor, material, financial resources, and time, and it is also easily affected by seasonality [20,21,22]. Although analytical methods can evaluate cutting forces and soil behavior such as soil displacement and velocity, the solutions of the analytical models are primarily based on the concept of equilibrium limits under quasi-static conditions and predetermined assumptions regarding the damaged surface, and the shape of the tillage tool must be simplified or of regular geometry [23,24]. The numerical simulation method has been developed as an efficient calculation method in recent years. Using numerical simulation methods to study the soil cutting process saves time and resources. It eliminates the interference of chance factors in the experimental research process and thus represents the modern design development direction [25,26,27].

The finite element method (FEM) [28,29,30], the discrete element method (DEM) [19,31,32,33], and computational fluid dynamics (CFD) [34,35,36] are the main numerical simulation methods used by scholars in the study of soil–tool interaction. In the classical finite element method (FEM), the nodes are fixed in the material, and the elements deform with the material. Large soil deformation, however, frequently results in severe mesh distortion, which significantly impacts the accuracy of the simulation results and can even cause the calculation to stop [37,38]. DEM is a good tool for simulating particle materials such as soil, but the stress–strain concept does not exist directly and cannot be calculated directly. DEM simulations require a rigorous and sophisticated method to calibrate the parameters representing soil mechanical properties. Furthermore, the large computational volume of the discrete element method (DEM) makes the computation very expensive [39,40]. Because the computational domain is fixed in space, the Eulerian description is used to prevent element deformation in computational fluid dynamics (CFD). Because the mesh is fixed in space, it is inapplicable when the material boundary is deformable. CFD simulation of the unconfined deformation of soil-free surfaces during cutting is therefore difficult [41].

The recently developed SPH method, also known as the smooth particle hydrodynamics method, is a meshless pure Lagrangian method [42]. Originally developed to solve astrophysical problems, the SPH method has since been extended to solve fluid flow, elastic–plastic flow, and brittle body damage [43,44]. In recent years, the SPH method has also been used for soil cutting. Man Hu et al. [45] investigated the cutting process and soil–tool interaction in non-cohesive and cohesive soils using an SPH model of soil–tool interaction based on an elastic–plastic structure, and the predicted cutting forces during tillage were in good agreement with the experimental results. As a result, using the SPH method to study soil–tool interactions is both practical and informative.

Few people have yet to apply the SPH method to simulate mouldboard plough–soil interactions. Thus, the following are the aims and objectives of the current work: (1) The mouldboard plough–soil interaction is simulated using the SPH method utilizing the finite element software LS-DYNA. (2) A coupled FEM-SPH model is developed to increase computational efficiency. (3) The simulation results are compared with those of the experiments [46] to verify the model’s reliability. (4) The effects of cutting speed, depth, inclination angle, and lifting angle on cutting resistance are investigated using the validated simulation model.

2. Materials and Methods

2.1. Brief Introduction to SPH

The SPH method is a meshless pure Lagrangian method with interpolation theory at its core [47,48], and its essence is to discretize a continuum with physical quantities such as mass and velocity into interacting particles and then to solve for the discrete individual particles to obtain the overall information of the continuum [49,50]. Any macroscopic variable function of a continuous variable field is represented as an integral form by smooth approximate interpolation, which can be expressed as:

$$u\left(x\right)={\displaystyle {\int }_{\Omega }u\left(\overline{x}\right)}\delta \left(x-\overline{x}\right)d\overline{x}$$

(1)

in which $u\left(x\right)$ is a continuous function, such as displacement, stress and so on, and $\delta \left(x-\overline{x}\right)$ is a delta function, which can be expressed as:

$$\delta \left(x-\overline{x}\right)=\left\{\begin{array}{l}1 x=\overline{x}\\ 0 x\ne \overline{x}\end{array}\right.$$

(2)

If a strongly peaked function $W\left(x-\overline{x},h\right)$ is employed to replace the delta function $\delta \left(x-\overline{x}\right)$, then Equation (1) can be rewritten as:

$$\langle u\left(x\right)\rangle ={\displaystyle {\int }_{\Omega }u\left(\overline{x}\right)}W\left(x-\overline{x},h\right)d\overline{x}$$

(3)

in which the strongly peaked function $W\left(x-\overline{x},h\right)$ is called the smooth kernel function; when $h\to 0$, the strongly peaked function tends to $\delta \left(x-\overline{x}\right)$. The most commonly used smooth kernel in SPH is the cubic B-sample, defined as:

$$\theta (\mu )=\left\{\begin{array}{l}(1-1.5{\mu }^2+0.75{\mu }^3)C (\left|\mu \right|\le 1)\\ 0.25C(2-\mu {)}^3 (1<\left|m\right|<2)\\ 0 (2\le \left|\mu \right|)\end{array}\right.$$

(4)

2.2. FEM–SPH Coupling Method

The methods in LS-DYNA that can be used to treat sliding and impact along the interface include the kinematic constraint method, the penalty method, the distributed parameter method, and others. As shown in Figure 1a, the SPH particle and the finite element in the coupling model established in this paper are divided into two parts: attachment coupling and contact coupling, with attachment coupling achieved using the kinematic constraint algorithm and contact coupling achieved using the penalty algorithm. All of the preceding algorithms are based on a master–slave scheme, in which the interfaces are defined as master and slave surfaces, and the nodes that lie within these surfaces are defined as master and slave nodes, respectively [51]. In general, it finds the nearest point in the master surface to each slave node to determine the penetration, unless the point is at the intersection of two master segments. As shown in Figure 1b, it is assumed that a master segment has been found for the slave node $n_s$ and that $n_s$ has not been determined as being located at the intersection of two master segments. The contact point is then defined as the nearest point on the master segment to the slave node $n_s$.

A kinematic constraint algorithm is used to achieve the attachment between the SPH particles and the finite element, which takes place at their interface [51]. During the computation, the slave nodes (SPH particles) are constrained to move on the master surface (finite element) after the impact, and the orthogonal projection keeps the relative position of the interpolated contact points on the master segment constant at each step of the computation. The tied interfaces are independently checked and updated at each time step. First, the nodal forces and masses of each slave node are distributed to the master node, and the mass increment and force increment of the master node define the segment containing the contact points, which can be calculated using Equations (5) and (6):

$$\Delta M_m^i={\varphi }_i\left({\xi }_c,{\eta }_c\right)M_s$$

(5)

$$\Delta f_m^i={\varphi }_i\left({\xi }_c,{\eta }_c\right)f_s$$

(6)

where ${\varphi }_i\left({\xi }_c,{\eta }_c\right)$ is the interpolation function, and $i$ is the number of master nodes; the force and mass increments calculated by Equations (5) and (6) will then be added to the force and mass vectors of the master surface. After completing the summation for all slave nodes, the acceleration of the $j$th node of the master surface $a_i^j$ can be calculated. Using $a_i^j$, Equation (7) can be used to calculate the acceleration $a_{is}$ of each slave node, updating the slave nodes’ velocity and displacement.

$$a_{is}={\displaystyle \sum _{j=1}^4{\varphi }_j}\left({\xi }_c,{\eta }_c\right)a_i^j$$

(7)

The penalty method is used to achieve contact between the finite element plough and the SPH area of the soil [51]. During the calculation, the contact will only occur if the penetration is positive; so, each slave node $n_s$ will be checked for penetration $l$ through the master surface.

$$l=n\times \left[g-r\left({\xi }_c,{\eta }_c\right)\right]<0$$

(8)

where $n$ is the vector normal to the master segment at contact point $\left({\xi }_c,{\eta }_c\right)$. $g$ and $r$ are the position vectors drawn to the slave and master nodes, respectively.

If penetration $l$ occurs $\left(l<0\right)$, an interface force vector $f_s$, which is proportional to the magnitude of the penetration, is applied between the slave node and its corresponding contact point.

$$f_s=-l{k}_i{n}_i$$

(9)

Then, the interface force applied to the four nodes $\left(i=1,2,3,4\right)$ of the master segment is expressed as:

$$f_m^i={\varphi }_i\left({\xi }_c,{\eta }_c\right)f_s$$

(10)

where ${\varphi }_i\left({\xi }_c,{\eta }_c\right)$ is the interpolation function.

$$k_i=\frac{f_{si}{K}_i{A}_i^2}{V_i}$$

(11)

where $k_i$ is the stiffness factor for master segment, $f_{si}$ is the scale factor, $K_i$ is the bulk modulus, $A_i$ is the surface area of the element, and $V_i$ is the volume.

2.3. Soil Constitutive Law

Over 300 material models are available in the LS-DYNA material library, but only a few of them have proven to be effective in describing flow and sliding properties. Remarkably, the elastic–plastic material MAT_010 (MAT_ELASTIC_PLASTIC_HYDRO) is suitable for simulating a wide range of materials, including fluid materials with large deformations [52]. As a result, in this paper, the elastic–plastic material MAT_010 is used to simulate soil. The type of soil is clay loam, and the parameters for the plough soil cutting simulation are shown in

Table 1. Simulation parameters of soil–plough coupling system [46,53].

Parameters	Value	Parameters	Value
Density of soil$/\left(\mathrm{k}\mathrm{g}/{\mathrm{m}}^3\right)$	2600	Static friction coefficient between soil and plough	0.6
Shear modulus of soil$/\mathrm{P}\mathrm{a}$	$1.2\times {10}^6$	Dynamic friction coefficient between soil and plough	0.1
Poisson’s ratio of soil	0.38	Static friction coefficient between two soil particles	0.422
Density of plough$/\left(\mathrm{k}\mathrm{g}/{\mathrm{m}}^3\right)$	7970	Dynamic friction coefficient between two soil particles	0.282
Young’s modulus of plough$/\mathrm{P}\mathrm{a}$	$2.10\times {10}^{11}$	Radius of particle$/\left(\mathrm{m}\mathrm{m}\right)$	5
Poisson’s ratio of plough	0.32	Gravity$/\left(\mathrm{m}/{\mathrm{s}}^2\right)$	9.806

.

If the effective stress (ES) and effective plastic strain (EPS) are not defined in MAT_010, the yield stress and plastic hardening modulus are taken from SIG0 and EH. In this case, the bilinear stress–strain curve shown in Figure 2 is obtained using hardening parameters. The yield strength is calculated as:

$${\sigma }_y={\sigma }_0+E_h{\overline{\epsilon }}^p+(a_1+p{a}_2)\max \left[p,0\right]$$

(12)

where $p$ is the positive compression pressure, and $E_h$ is the plastic hardening modulus defined by Young’s modulus $E$ and the tangent modulus $E_t$, which are defined as follows:

$$E_h=\frac{E_{t}E}{E-E_t}$$

(13)

If the effective stress (ES) and effective plastic strain (EPS) are specified, a curve similar to the one shown in Figure 2 can be defined. The effective stress is defined by the deviatoric stress tensor $s_{ij}$ as:

$$\overline{\sigma }={\left(\frac{3}{2}{s}_{ij}{s}_{ij}\right)}^{{}^1{/}_2}$$

(14)

and effective plastic strain by:

$${\overline{\epsilon }}^p={{\displaystyle {\int }_0^t\left(\frac{2}{3}{D}_{ij}^p{D}_{ij}^p\right)}}^{{}^1{/}_2}dt$$

(15)

where $t$ is the time, and $D_{ij}^p$ is the plastic component of the deformation tensor rate. In this case, ignoring the plastic hardening modulus, the yield stress is given as:

$${\sigma }_y=f\left({\overline{\epsilon }}^p\right)$$

(16)

where the value of $f\left({\overline{\epsilon }}^p\right)$ is obtained by interpolation of the data curve.

3. Numerical Simulation

3.1. Establishment of the Finite Element Model of the Plough

In this paper, the plough chosen for the tractor CF700 tillage system in the experiment of Jinming Zhang et al. [46] is used as the research object. The numerical simulation method is used to study the entire process of soil cutting by the plough, and the numerical simulation results are compared with those of the experiment for analysis. This paper only studies the working condition of a single plough in the soil to simplify the model and improve calculation efficiency.

The plough is modeled utilizing 3D modeling software using the horizontal straight element line method, and its geometric model is shown in Figure 3a. Considering the irregularity of the plough’s surface, the geometric model of the plough is meshed, and the finite element mesh is free-meshed with a mesh size of 10 mm. According to the literature [46,53], the plough chosen for this study is made of 65 Mn steel with a material density of $7970 \mathrm{k}\mathrm{g}/{\mathrm{m}}^3$, a modulus of elasticity of $2.10\times {10}^{11} \mathrm{P}\mathrm{a}$, and a Poisson’s ratio of 0.32. The model was imported into LS-DYNA; the element type was set to Solid; the material model was set to MAT_RIGID; and the constraints were set to move in the y and z directions and rotate in the x, y, and z directions. Figure 3b depicts the plough’s finite element model.

3.2. Establishment of SPH Model, FEM–SPH Coupling Model

3.2.1. SPH Model

There are two methods in LS-DYNA for generating SPH soil particles: one is to convert the established finite element mesh into SPH particles; this method is mostly used for modeling objects with irregular and complex structures. The other is to generate SPH particles directly; this is the best way to model objects with simple structures. The direct SPH particle generation method is used in this paper for soil modeling. First, the established finite element mesh model of the plough was imported into the LS-PrePost. The rectangular soil model was generated after setting the particle generation approach, the number of particles in each direction, the filling rate, the particle density, and other parameters in the Mesh_Sph Generation, with a soil size of 3 × 2 × 0.6 m and particle number of 132,000. Second, parameters such as particle type and radius were set in Settings_General Settings. Finally, Element Tools_Transform was used to adjust the cutting position of the plough based on the plough’s cutting mode and the boundary condition processing requirements. Figure 4 depicts the SPH model of the plough soil cutting.

3.2.2. FEM-SPH Coupling Model

The soil in the coupled model comprises two parts: the finite element mesh and the SPH particles, with the SPH particles obtained through finite element mesh transformation. First, the established finite element mesh model of the plough is imported into the LS-PrePost, and Mesh_ShapeM generates a finite element mesh with a soil size of 3 × 2 × 0.6 m and an element size of 0.024 m. Second, by using Mesh_Sph Generation, the finite element mesh of the middle part of the soil is transformed into SPH particles with a soil size of 3 × 0.96 × 0.46 m. Finally, in Settings_General Settings, parameters such as particle type and particle radius are set, and the cutting position of the plough is adjusted using Element Tools_Transform. Figure 5 depicts the FEM–SPH plough soil cutting model.

3.2.3. Boundary Conditions Imposed

The plough’s depth into the soil during the soil cutting process is 0.1 m, and it moves forward with the tractor. To reduce invalid cutting time, the plough should be as close to the soil edge as possible. In the numerical simulation, the following boundary conditions are imposed: ① Select the numerical algorithm as the SPH algorithm. ② The plough’s forward speed is 1 m/s, and the cutting depth is 0.1 m. ③ The soil material was chosen as MAT_010 (MAT_ELASTIC_PLASTIC_HYDRO). ④ Using SPC, set the soil bottom, with the left and right sides as fixed constraints, and apply gravity to the plough and the soil. ⑤ Set the contact between the plough and the soil as AUTOMATIC_NODES_TO_SURFACE. Set the contact between FEM and SPH soil in the coupled model as a TIED contact (TIED_NODES_TO_SURFACE). ⑥ Set the solution time to 4 s. Complete the other keyword settings in the LS-PrePost and generate the K file to be submitted to the LS-DYNA solver for solution.

3.3. Parameter Study

To investigate the effect of cutting speed and depth on cutting resistance, the cutting resistance corresponding to the different cutting speeds (1 m/s, 2 m/s, 3 m/s, 4 m/s, and 5 m/s) under the different cutting depth conditions (0.1 m, 0.15 m, 0.2 m, 0.25 m, and 0.3 m) of the plough was simulated using the simulation model.

The inclination angle $\alpha$ and lifting angle $\beta$ of the plough are depicted in Figure 6 and are as given in the study in [54]. The cutting resistance at different cutting depths was simulated using a simulation model to investigate the effect of plough inclination and lifting angles on cutting resistance. For the inclination angle, the cutting resistance was simulated for different depths (0.15 m, 0.2 m, and 0.25 m) as the inclination angle was increased from $30°$ to $75°$. For the lifting angle, the cutting resistance was simulated for different depths (0.15 m, 0.2 m, and 0.25 m) as the lifting angle was increased from $15°$ to $60°$. In the simulation, the plough’s cutting speed was 1 m/s.

4. Results and Discussion

4.1. Cutting Process and Cutting Resistance Analysis

Figure 7 depicts the soil cutting process and the cutting resistance using the SPH model and the coupled FEM-SPH model, which clearly shows the trajectory and displacement changes of the soil during the cutting process. When the plough moves into contact with the soil, the soil is destroyed due to shearing and squeezing (Figure 7a,d); as the plough cuts entirely into the soil, the contact area between the plough and the soil increases, and the soil is gradually destroyed to a greater extent (Figure 7b,e). When the plough moves forward to the edge of the soil and then leaves the soil, the soil remains in motion due to inertia, causing the plough to cut the soil continuously (Figure 7c,f). The soil cutting processes in both simulation models are very stable, and both can reflect the soil being turned up and falling. The soil cutting process is consistent with the actual situation, indirectly confirming the simulation model’s reliability.

The cutting resistance of a plough is the reaction force of the soil on the plough during the cutting process. Loading the rcforc file into the LS-DYNA displays the plough’s cutting resistance. Figure 7 shows the combined cutting resistance curves obtained from the simulations of the SPH model and the coupled FEM-SPH model. As shown in the figure, there is no initial contact between the plough and the soil, and the cutting resistance value is 0. The cutting resistance increases gradually as the plough contacts and cuts into the soil. At 0.8 s, the plough completely cuts into the soil, and the cutting resistance reaches and fluctuates near its maximum value. The plough starts cutting out of the soil at 3 s, and the cutting resistance gradually decreases. At 4 s, the plough completely removes the soil, and the cutting resistance value is 0.

The plough cut completely into the soil at 0.8 s and began to cut out of the soil at 3 s. The cutting resistance data from the two numerical simulation models were extracted from 0.8 s to 3 s. After analysis and calculation, the average cutting resistances of the SPH and coupled FEM-SPH models were 349.48 N and 306.25 N, respectively.

4.2. Stress Variation of Soil

Figure 8 depicts the change in soil effective stress during the plough soil cutting process using the coupled FEM-SPH model. The effective stress on the soil increases as the contact area between the plough and the soil increases during the soil cutting process. When the soil is completely destroyed, the effective stress reaches a stable value. During the soil cutting process, the maximum effective stress on the soil is 4000 Pa. As can be seen from the entire cutting process, the soil is cut relatively steadily, and the effective stress on the soil varies less.

4.3. Model Validation

To validate the reliability of the simulation model, the simulation results of the SPH model and the coupled FEM-SPH model were compared with those of the experiment and the DEM simulation [46]. Figure 9 depicts the simulation and experimental cutting resistance curves at 1 m/s cutting speed and 0.1 m depth of cut (Figure 9a shows the horizontal cutting resistance of the plough, and Figure 9b shows the vertical cutting resistance of the plough.). Figure 9 shows that the cutting resistance of the two numerical simulation models is basically consistent with that of the experiment, proving the model’s reliability.

When the plough makes contact with the soil, the cutting resistance in the horizontal and vertical directions increases; it stabilizes after a while and fluctuates up and down near the stable value. The following can be derived from both the simulation model and the literature: The SPH model’s average cutting resistances in the horizontal and vertical directions are 124.97 N and 62.50 N, respectively. The coupled FEM-SPH model’s average cutting resistances in the horizontal and vertical directions are 130.07 N and 72.62 N, respectively. The experiment’s average cutting resistances in the horizontal and vertical directions are 122.32 N and 64.87 N, respectively. The DEM model’s average cutting resistances in the horizontal and vertical directions are 121.05 N and 62.64 N, respectively. They are calculated as follows: The relative errors between the SPH model and the experiment are 2.17% and 3.65%, respectively. The relative errors between the coupled FEM–SPH model and the experiment are 5.96% and 10.67%, respectively. The relative errors between the DEM model and the experiment are 1.04% and 3.43%, respectively.

Regarding computational efficiency and accuracy, the SPH model and the coupled FEM-SPH model have advantages and disadvantages. Regarding computational efficiency, the SPH model takes 5 h and 2 min to compute using the same computing device, while the FEM-SPH coupled model takes 38 min. The FEM-SPH coupled model is more computationally efficient than the SPH model. In terms of computational accuracy, it is obvious from the above error analysis that the SPH model has a higher computational accuracy than the coupled FEM-SPH model. The DEM model has the lowest error and the highest accuracy. However, the large computational volume of DEM may render it inefficient in terms of computational efficiency.

Figure 10 depicts the linear regression plots of cutting resistance in the horizontal and vertical directions for both the experiment and the simulation. As shown in the figure, both simulation models predicted a cutting resistance with good regression results. Numerically, the SPH model has the best regression results for cutting resistance in the horizontal direction, followed by the FEM-SPH model in the vertical direction, then the SPH model in the vertical direction and, finally, the FEM-SPH model in the horizontal direction.

4.4. Effect of Cutting Speed and Depth on Cutting Resistance

Figure 11 depicts the variation in cutting resistance as the plough’s cutting speed increases from 1 m/s to 5 m/s under various cutting depth (0.1 m, 0.15 m, 0.2 m, 0.25 m, and 0.3 m) operating conditions (Figure 11a depicts the horizontal cutting resistance; Figure 11b depicts the vertical cutting resistance; and Figure 11c depicts the combined cutting resistance). Fitting the data into the graph results in the following: The cutting resistance increases with the increase in cutting speed and has a quadratic relationship under different cutting depth conditions. The cutting resistance increases with the depth of the cut and has a quadratic relationship under different cutting speeds. According to the study in [55], the relationship between cutting resistance and cutting speed and depth can be expressed as:

$$f=a{v}^2+b$$

(17)

$$f=p{d}^2+qd$$

(18)

where $a$ and $b$ are empirical functions of the depth $d$, and $p$ and $q$ are empirical functions of the velocity $v$, which depends on specific soil properties and geometric parameters of the plough. Obviously, cutting resistance is a quadratic function of cutting speed and depth. Therefore, the results predicted by the simulations in this study are consistent with those of the above studies.

4.5. Effect of Inclination Angle $\alpha$ and Lifting Angle $\beta$ of the Plough on Cutting Resistance

Figure 12 depicts the variation in cutting resistance as the inclination angle increases from $30°$ to $75°$ at a cutting speed of 1 m/s and different depths of cut (0.15 m, 0.2 m, and 0.25 m) (Figure 12a depicts the horizontal cutting resistance; Figure 12b depicts the vertical cutting resistance; and Figure 12c depicts the combined cutting resistance). It can be seen that under different cutting depth operating conditions, the horizontal cutting resistance and combined cutting resistance increase with the increase in inclination angle $\alpha$, whereas the vertical cutting resistance remains essentially constant with the increase in inclination angle $\alpha$.

Figure 13 depicts the variation in cutting resistance as the lifting angle increases from $15°$ to $60°$ at a cutting speed of 1 m/s and different depths of cut (0.15 m, 0.2 m, and 0.25 m) (Figure 13a depicts the horizontal cutting resistance; Figure 13b depicts the vertical cutting resistance; and Figure 13c depicts the combined cutting resistance). It can be seen that, under different cutting depth operating conditions, horizontal cutting resistance and combined cutting resistance increase with the increase in lifting angle $\beta$, whereas vertical cutting resistance decreases with the increase in lifting angle $\beta$.

5. Conclusions

The mouldboard plough-soil interaction was simulated using the SPH method and the FEM-SPH coupling method. The cutting resistance obtained from the two simulation models was compared with that of the experiment in the literature, and the curves of the simulation and experiment were in good agreement, verifying the reliability of the SPH model and the FEM-SPH coupled model. As a result, the two simulation models can be used as a quick and accurate tool to predict cutting resistance during soil cutting.

The SPH model and the FEM-SPH coupled model can both accurately simulate the soil cutting process, but each has advantages and disadvantages in terms of computational efficiency and accuracy. The SPH model takes 5 h and 2 min to compute, while the FEM-SPH coupled model takes 38 min; obviously, the computational efficiency of the FEM-SPH coupled model is higher. The relative errors between the SPH model and the experiment are 2.17% and 3.65%, respectively. The relative errors between the FEM-SPH coupled model and the experiment are 5.96% and 10.67%, respectively. Obviously, the SPH model has a higher computational accuracy.

The maximum effective stress on the soil during the entire cutting process was 4000 Pa, and the maximum effective stress was primarily concentrated on the soil in contact with the mouldboard plough. The average cutting resistances predicted by the SPH model and the FEM-SPH coupled model, respectively, were 349.48 N and 306.25 N; these resistances are useful for precision tillage.

The validated simulation model was used to investigate the effects of cutting speed, depth of cut, inclination angle, and lifting angle of the mouldboard plough on cutting resistance. Simulation studies have shown that the cutting resistance increases with the increase in cutting speed and is quadratic; the cutting resistance increases with the increase in cutting depth and is quadratic; the horizontal cutting resistance and the combined cutting resistance increase with the increase in inclination angle, while the vertical cutting resistance remains essentially constant with the increase in inclination angle; the horizontal cutting resistance and combined cutting resistance increase as the lifting angle increases, while the vertical cutting resistance decreases as the lifting angle increases.