Simulation of Plastic Deformation Failure of Tillage Tools Based on the Smoothed Particle Hydrodynamics Method

Abstract

The problems of large deformations, failures, and fractures that agricultural tillage tools may encounter during the cultivation process has long been a concern in the field of agricultural machinery design and manufacturing. It is important to establish a more accurate numerical model to effectively predict tools’ plastic deformation failures and ductile fracture failures. This research develops a numerical model for predicting the plastic deformation failure and ductile fracture failure of agricultural tillage tools using the smoothed particle hydrodynamics (SPH) method and the Johnson–Cook constitutive model. The model uses the Drucker–Prager criterion to describe the elastic–plastic constitutive behavior of the soil, the von Mises criterion to describe the Johnson–Cook constitutive model of the tool, and the coupling condition with the Lennard-Jones repulsive force to describe the interaction between the tool and soil. The numerical results show that the proposed model can effectively simulate the interaction between the tool and soil, as well as the tool’s plastic deformation failure and ductile fracture failure during the agricultural cultivation process. It can also predict the variation trend of the cutting force of the tool. This helps to provide a new approach for the numerical simulation of such problems.

1. Introduction

Agricultural tillage tools (such as plowshares, subsoilers, disc harrows, rotary blades, etc.) [1] may encounter plastic deformation failure and ductile fracture failure during cultivation operations, especially in regions with poor soil conditions (such as the rocky and sticky soil in mountainous areas, where the probability of tool failure is high) [2]. Assessing and predicting the likelihood of plastic deformation failure and ductile fracture failure in agricultural tillage tools and evaluating the tool–soil interaction during agricultural tillage operations has always been a challenging issue in agricultural tillage tool research. This is because the condition of the tillage tool directly affects the reliability and efficiency of the overall tillage machinery operation as a key component [3,4]. Therefore, establishing a reliable mechanical model to study the cutting process and the potential large deformation failures of the tool would be beneficial. Such a model could assist in optimizing tool design and reducing the failure rate of tillage machinery.

Currently, research on agricultural tillage tools primarily includes two approaches: field experiments and numerical simulations. Field experiments are the major approach used for designing and studying agricultural tillage tools [5]. However, field experiments are limited to studying the cutting process and tool conditions using external macroscopic morphology and measurements. This makes it difficult to reflect the dynamic mechanism arising from the internal structure of the tool and soil through field experiments alone. Additionally, conducting comprehensive field experiments can be time-consuming and laborious [6]. With the rapid development of computer technology and numerical methods, simulating the dynamic responses of the tool and soil through numerical methods has become an effective approach to studying tool–soil interactions, which facilitates research into tool optimization, failure prediction, and tillage strategy [4].

In the study of tillage tool–soil interactions in the agricultural field, tools are typically simplified and modeled as either rigid or deformable bodies. Currently, the assumption of tool rigidity is more widely used. This simplification does not consider tool deformation, greatly reducing the modeling and computational complexity [7,8,9,10,11,12,13]. For example, Yang et al. [14] simulated the soil-cutting process using the finite element method (FEM) and studied the cutting force under different parameters, such as tool rake angles, tillage depths, and friction coefficients; however, they did not model tool deformation. Similarly, Fang et al. [15] verified a discrete element method (DEM) tool–soil interaction model by establishing an experimental setup and analyzing the changes in the resultant force, horizontal force, and vertical force of the rotary blade; this also assumed the tool to be a rigid body. Gao et al. [7] developed a soil–tool interaction model based on SPH, verified the accuracy of the soil–tool interaction model; they also assumed a rigid tool. In recent years, some researchers have started to model tools as deformable bodies using multi-physics joint simulations, obtaining information such as tool stress and deformation. However, this approach requires significant computational resources compared with the assumption of rigidity. Zhang et al. [16] conducted FEM studies on the maximum stress value, maximum total deformation, and minimum fatigue life of a rotary blade. The simulation results agreed with the actual conditions, verifying the feasibility of FEM simulations for small deformations of tools. Zhang et al. [17], focusing on the easy wear failure problem of soil-engaging components in agricultural machinery, used DEM and FEM for a combined simulation. They obtained parameters such as the average deformation and equivalent stress after a sample wear test, further evaluating the wear characteristics and exploring the wear-resistance mechanism.

In general, the aforementioned methods can simulate the cultivation process and analyze tool forces to a certain extent. However, simulating large or fracturing deformations under certain scenarios remains difficult, especially when trying to predict failure by considering tools as deformable bodies with material properties. The classic FEM is prone to numerical divergence caused by mesh distortion when processing large deformations [4]. DEM can simulate large deformation problems, but it is more suitable for materials with non-cohesive soil particles, and the calibration of soil parameters is complicated [18]. Some coupled methods require continuous updating of the grid during the calculation process. As a meshless method, the SPH method, with its superior adaptivity, allows the solution process to remain unaffected by the randomness of particle distribution; this can handle extremely large deformation problems well. It also has a wide range of applications in solving highly nonlinear problems such as high-speed collisions and explosions [13]. Therefore, the SPH method has good potential for modeling deformable tool–soil interactions.

This paper establishes a coupled numerical model based on the smoothed particle hydrodynamics (SPH) method, to consider the interaction between deformable tools and soil. The model utilizes the Drucker–Prager constitutive relation to describe the elastoplastic deformation behavior of the soil. It also implements the Johnson–Cook constitutive model to characterize the progressive plastic deformation and fracture evolution of the tool under loading. Through simulation and analysis, the model is found to accurately predict the nonlinear tool–soil interaction as well as the complex plastic deformation and fracture processes in the tool. Compared with traditional methods assuming rigid tools, this coupled model can more accurately describe the nonlinear tool–soil interaction and complex deformation mechanisms that occur. By accounting for large deformations and potential failure modes, this study provides a new perspective for the optimization design of agricultural equipment and tools.

2. Simulation Approaches

2.1. The Basics of SPH Method

The basic idea of SPH can be summarized as follows: use a group of particles with physical properties to discretize the problem domain, approximate the governing equations of the fluid to the form of the ordinary differential SPH equations, and solve these using the weighted averages of the particles in the support domain to obtain the desired physical quantities. There are two key steps in the construction of the SPH equation: kernel function approximation and particle approximation.

2.1.1. Kernel Function Approximation

The integral expression of the function f(x) at point x in the problem domain $\Omega$ is defined as

$$f\left(x\right)={{\displaystyle \int }}_{\Omega }^{}f\left(x\prime \right)\delta \left(x-x\prime \right)dx\prime$$

(1)

where $\delta \left(x-x\prime \right)$ is the Dirac function:

$$\delta \left(x-x\prime \right)=f\left(x\right)=\{\begin{array}{c}\infty , x=x\prime \\ 0, x\ne x\prime \end{array}$$

(2)

By replacing the Dirac function with a smoothing kernel function, the expression of the h kernel function can be approximated as

$$\langle f\left(x\right)\rangle ={{\displaystyle \int }}_{\Omega }^{}f\left(x\prime \right)W\left(x-x^{\prime },h\right)dx\prime$$

(3)

In this equation, $W\left(x-x^{\prime },h\right)$ represents the smoothing kernel function; h is the smooth length, which is used to define the influence of the smoothing kernel function [19]. This study selects the cubic spline smoothing function [20] as the kernel function for SPH research:

$$W\left(R,h\right)={\alpha }_d\times \{\begin{array}{c}\frac{2}{3}-R^2+\frac{1}{2}{R}^3, 0\le R<1\\ \frac{1}{6}{\left(2-R\right)}^3, 1\le R<2\\ 0, R\ge 2\end{array}$$

(4)

In this equation, $R=\left|x - x’\right|/h$, for two-dimensional spaces, and ${\alpha }_d$ takes values of 15/(7$\pi$h²).

2.1.2. Particle Approximation

The particle approximation discretizes the problem domain into a finite number of particles with independent spaces. The field variable on each particle can be obtained by summing the weighted field variables of all neighboring particles within the support domain [21]. In the SPH method, the volume of particle j is represented by using the volume ${\Delta V}_j$ to replace the infinitesimal element dx in Equation (3); thus, it can be expressed as [19]

$$\Delta V_j=m_j/{\rho }_j$$

(5)

In this equation, $m_j$ is the mass of particle j (j = 1, 2, …, N), where N is the total number of particles within the support domain of particle j; ${\rho }_j$ is the density of particle j [19]. Substituting Equation (5) into Equation (3) yields an approximate equation for particle i:

$$\langle f\left(x_i\right)\rangle ={\displaystyle \sum }_{j=1}^N\frac{m_j}{{\rho }_j}f\left(x_j\right)·W\left(x_i-x_j,h\right)$$

(6)

In this equation, $x_i$ and $x_j$ represent the respective position function variables of particle i and particle j.

2.2. Discretization of Governing Equations

The mass conservation equation in the Navier–Stokes equations under the SPH method is

$${\rho }_i=\frac{{{\displaystyle \sum }}_{j=1}^N{m}_j{W}_{ij}}{{{\displaystyle \sum }}_{j=1}^N\frac{m_j}{{\rho }_j}{W}_{ij}}$$

(7)

The momentum conservation equation in the Navier–Stokes equations under the SPH method is

$$\frac{D{v}_i^{\alpha }}{Dt}={\displaystyle \sum }_{j=1}^N{m}_j\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}\right)\frac{\partial W_{ij}}{\partial x_i^{\beta }}+f^{\alpha }$$

(8)

The energy conservation equation in the Navier–Stokes equations under the SPH method is

$$\frac{D{e}_i}{Dt}=\frac{1}{2}{\displaystyle \sum }_{j=1}^N{m}_j\left(\frac{P_i+P_j}{{\rho }_i{\rho }_j}\right)v_{ij}^{\beta }\frac{\partial W_{ij}}{\partial x_i^{\beta }}+\frac{{\mu }_i}{2{\rho }_i}{\epsilon }_i^{\alpha \beta }{\epsilon }_j^{\alpha \beta }$$

(9)

In the equation, v represents velocity, e represents internal energy, $\sigma$ represents the stress tensor, P represents pressure, m represents mass, $\rho$ represents density, $\mu$ represents the dynamic viscosity coefficient, $\epsilon$ represents the strain rate, N represents the total number of particles within the support domain, x represents the particle coordinate position, and $\alpha \mathrm{and} \beta$ denote Cartesian components; i and j represent different particles.

In order to avoid the issues of stress instability and nonphysical oscillations in traditional SPH methods [22], the widely used artificial stress term $K_{ij}$ and artificial viscosity term ${\Pi }_{ij}$ [23] are employed in the solution. The momentum from Equation (8) is modified as

$$\frac{D{v}_i^{\alpha }}{Dt}={\displaystyle \sum }_{j=1}^N{m}_j\left(\frac{{\sigma }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\sigma }_j^{\alpha \beta }}{{\rho }_j^2}+K_{ij}+{\Pi }_{ij}\right)\frac{\partial W_{ij}}{\partial x_i^{\beta }}+f^{\alpha }$$

(10)

2.3. Constitutive Model of the Soil

The soil particle stress rate equation in the SPH form is constructed based on an elastoplastic constitutive model adopting the Drucker–Prager (D–P) yield criterion for soil particles:

$$\frac{D{\sigma }_i^{\alpha \beta }}{Dt}={\sigma }_i^{\alpha \gamma }{\dot{\omega }}_i^{\beta \gamma }+{\sigma }_i^{\gamma \beta }{\dot{\omega }}_i^{\alpha \gamma }+2G{\dot{e}}_i^{\alpha \beta }+K{\dot{\epsilon }}_i^{\gamma \gamma }{\dot{\delta }}_i^{\alpha \beta }-{\dot{\lambda }}_i\left[9K\sin \psi {\delta }^{\alpha \beta }+\frac{G}{\sqrt{J_2}}{s}_i^{\alpha \beta }\right]$$

(11)

In this equation, the first two terms are the Jaumann stress rate [23], ${\dot{\epsilon }}^{\alpha \beta }$ represents the strain rate tensor, ${\dot{\omega }}^{\alpha \beta }$ represents the rotation rate tensor, and ${\dot{\lambda }}_i$ represents the plastic multiplier; these are defined as

$${\dot{\epsilon }}_i^{\alpha \beta }=\frac{1}{2}\left[{\displaystyle \sum }_{j=1}^N\frac{m_j}{{\rho }_j}\left(v_j^{\alpha }-v_i^{\alpha }\right)\frac{\partial W_{ij}}{\partial x_i^{\beta }}+{\displaystyle \sum }_{j=1}^N\frac{m_j}{{\rho }_j}\left(v_j^{\beta }-v_i^{\beta }\right)\frac{\partial W_{ij}}{\partial x_i^{\alpha }}\right]$$

(12)

$${\dot{\omega }}_i^{\alpha \beta }=\frac{1}{2}\left[{\displaystyle \sum }_{j=1}^N\frac{m_j}{{\rho }_j}\left(v_j^{\alpha }-v_i^{\alpha }\right)\frac{\partial W_{ij}}{\partial x_i^{\beta }}-{\displaystyle \sum }_{j=1}^N\frac{m_j}{{\rho }_j}\left(v_j^{\beta }-v_i^{\beta }\right)\frac{\partial W_{ij}}{\partial x_i^{\alpha }}\right]$$

(13)

$${\dot{\lambda }}_i=\frac{3{\alpha }_{\varphi }K{\dot{\epsilon }}_i^{\gamma \gamma }+\left(G/\sqrt{J_2}\right)s_i^{\alpha \beta }{\dot{\epsilon }}_i^{\alpha \beta }}{27{\alpha }_{\varphi }K\sin \psi +G}$$

(14)

2.4. Constitutive Model of the Cutting Tool

The total stress tensor ${\sigma }^{\alpha \beta }$ in Equation (10) is composed of the volumetric pressure P and the deviatoric shear stress ${\tau }^{\alpha \beta }$:

$${\sigma }^{\alpha \beta }=-P{\delta }^{\alpha \beta }+{\tau }^{\alpha \beta }$$

(15)

By substituting Equation (15) into Equation (10), we obtain

$$\frac{D{v}_i^{\alpha }}{Dt}={\displaystyle \sum }_{j=1}^N{m}_j\left[-\left(\frac{P_i}{{\rho }_i^2}+\frac{P_j}{{\rho }_j^2}\right){\delta }^{\alpha \beta }+\left(\frac{{\tau }_i^{\alpha \beta }}{{\rho }_i^2}+\frac{{\tau }_j^{\alpha \beta }}{{\rho }_j^2}\right)+K_{ij}+{\Pi }_{ij}\right]\frac{\partial W_{ij}}{\partial x_i^{\beta }}+f^{\alpha }$$

(16)

2.4.1. Equation of State

The volumetric pressure P of the tool in Equation (16) is calculated using the equation of state, as

$$P=\frac{{\rho }_0{C}_0^2\eta \left[1+\left(1-\frac{{\Gamma }_0}{2}\right)\eta \right]}{\left[1-\left(S_{a }-1\right)\eta \right]}+{\rho }_0{\Gamma }_{0}e$$

(17)

In this equation, $\eta =\frac{\rho }{{\rho }_0} - 1$, ${\Gamma }_0$ is the Gruneisen’s constant, $S_a$ is a linear Hugoniot slope coefficient, and $e$ is the internal energy per unit of mass.

2.4.2. von Mises Yield Criterion

The shear stress $\tau$ of elastic solids can be obtained from

$$\frac{\mathrm{d}{\tau }_i^{\alpha \beta }}{\mathrm{d}t}={\tau }_i^{\alpha \gamma }{\dot{\omega }}_i^{\beta \gamma }+{\tau }_i^{\gamma \beta }{\dot{\omega }}_i^{\alpha \gamma }+2G\left({\dot{\epsilon }}_i^{\alpha \beta }-\frac{1}{3}{\delta }^{\alpha \beta }{\dot{\epsilon }}_i^{\gamma \gamma }\right)$$

(18)

In this equation, G represents the shear modulus of the material, ${\dot{\epsilon }}^{\alpha \beta }$ denotes the strain rate tensor, and ${\dot{\omega }}^{\alpha \beta }$ represents the rotation rate tensor, as shown in Equations (12) and (13).

According to [24], based on the von Mises yield criterion, the aforementioned elastic constitutive relationship can be extended to plastic behavior:

$$f_y=\frac{{\sigma }_y}{\sqrt{3J_2}}$$

(19)

In this equation, ${\sigma }_y$ represents the yield stress, and $J_2$ denotes the second stress tensor invariant. When the von Mises yield criterion is satisfied ($f_y < 1$), the material is in the plastic state; then, the shear stress is brought back to the yield surface. The shear stress ${\tau }^{\alpha \beta }\prime$ at this point is

$${\tau }^{\alpha \beta }\prime =f_y{\tau }^{\alpha \beta }$$

(20)

2.4.3. Johnson–Cook Constitutive Model

The Johnson–Cook constitutive model can describe the plastic deformation behavior of metals under conditions of large deformations, high strain rates, and elevated temperatures [25]. In this study, the Johnson–Cook constitutive model is selected as the material model for the cutting tool, where the yield stress of this constitutive model is represented as

$${\sigma }_y=\left[A+B{\left({\epsilon }_{eff}^p\right)}^N\right]\left[1+Cln\left(\frac{\dot{{\epsilon }_{eff}^p}}{\dot{{\epsilon }_0}}\right)\right]\left[1-{\left(T^{*}\right)}^M\right]$$

(21)

In this equation, ${\epsilon }_{\mathrm{eff}}^p(={{\displaystyle \int }}^{ }\dot{{\epsilon }_{\mathrm{eff}}^p}dt)$ represents the equivalent plastic strain; ${\dot{\epsilon }}_{eff}^p\left(=\sqrt{\frac{2}{3}{\dot{\epsilon }}^{\alpha \beta }{\dot{\epsilon }}^{\alpha \beta }}\right)$ represents the equivalent plastic strain rate; $\dot{{\epsilon }_0}$ represents the reference value of the equivalent plastic strain rate; A, B, C, N, and M are constants related to the material; and $T^{*}$ represents the normalized temperature [24]:

$$T^{*}=\frac{T-T_{ref}}{T_{melt}-T_{ref}}$$

(22)

In this equation, $T_{ref}$ represents the reference temperature and $T_{melt}$ represents the melting temperature of the corresponding material.

The Johnson–Cook constitutive model can be combined with progressive damage and failure models to specify different damage initiation criteria and damage evolution laws, to demonstrate the phenomenon of fracture [25].

The progressive damage and failure criterion is employed to simulate the influence of tool deformation failure. Parameter D is used to measure the local damage state, as

$$D={{\displaystyle \sum }}^{ }\frac{\Delta {ϵ}_{eff}^p}{{ϵ}_{failure}}$$

(23)

In this equation, $\Delta {\epsilon }_{eff}^p$ represents the increment of equivalent plastic strain that occurs within an integration cycle, and ${\epsilon }_{failure}$ denotes the failure strain:

$${\epsilon }_{failure}=\left[D_1+D_2\exp (D_3{\sigma }^{*}\right]\left[1+D_4\ln \left(\frac{\dot{{\epsilon }_{eff}^p}}{\dot{{\epsilon }_0}}\right)\right]\left[1+D_5{T}^{*}\right]$$

(24)

When the parameter D in Equation (23) exceeds 1 for particle i, a material failure occurs and the corresponding shear stress at point i is set to zero [24].

2.5. Tool–Soil Interaction Model

This paper introduces a cutting tool–soil interaction model that employs the soil–solid interface coupling condition [26], which satisfies the following conditions for each particle located at the interface:

$$v_{soil}^{\alpha }=v_{solid}^{\alpha }$$

(25)

$${\sigma }_{soil}^{\alpha \beta }={\sigma }_{solid}^{\alpha \beta }$$

(26)

To satisfy the conditions above, particles near the cutting tool–soil interface should meet the following requirements: During the SPH particle weighted summation process, for all particles within a distance of 2 h from the interface, we should neglect the differences in particle types (i.e., soil particles and tool particles), as shown in Figure 1. In other words, the velocity tensors and stress tensors of the tool particles can be applied to the soil particles through a smoothing function, and vice versa [26].

To prevent mutual penetration between the soil and tool particles, a penalty force similar to the Lennard-Jones repulsive force [27] is applied to them when particles of different materials approach each other near the interface.

2.6. Model Implementation Process

Based on the above theory, this paper uses Fortran to write an SPH program for numerical computation of the tool-cutting process. The implementation process of the program is shown in Figure 2. The main task in establishing the initial model is to write the initial coordinates, density, velocity, and force conditions of soil particles and tool particles. Then, using the elastic–plastic constitutive model and the Johnson–Cook constitutive model, and combining them with a series of numerical processing methods, the internal force parameters of soil particles and tool particles are calculated separately, and the interaction is solved by coupling conditions. Finally, using temporal integration, the particle information is updated, and the desired data are outputted.

3. Numerical Example and Result Analysis

3.1. Model Construction and Parameter Setting

This paper establishes a two-dimensional numerical model to study the progressive deformation and fracturing behavior of tillage tools under different representative working conditions. As depicted in Figure 3a, the proposed model consists of a deformable tillage tool tilling a simplified rectangular soil domain. The tillage tool is modeled with dimensions of length 25 mm and thickness 1.5 mm, and it tills the soil at an angle of 70° to a depth of 10 mm. To enable clear observations of the dynamic cutting process, the tillage tool is modeled using a relatively soft aluminum alloy rather than the wear-resistant steels commonly utilized in actual tillage operations. While unrealistic in terms of material properties, the use of a lower-strength aluminum alloy serves to amplify the progressive deformation response and ductile fracture characteristics of the tillage tool, providing more prominent visualizations of the tool’s behavior. The soil domain is modeled as a simple fixed rectangular region with dimensions of length 50 mm and height 15 mm and is discretized into a network of uniformly spaced particles with a spacing of 1 mm between particles. This provides an appropriate balance between computational efficiency and solution accuracy. Typical parameters for a representative cohesive agricultural soil are utilized, accounting for plasticity and cohesion. The thermodynamic and damage parameters for the aluminum alloy are also selected to accurately characterize the material’s progressive ductile deformation and fracture behavior under the complex loading imparted by the soil-cutting process. As summarized in

Table 1. Material parameters.

Material Type	Material Properties	Symbol	Value
Soil	Density	$\rho$	1540 kg/m3
	Cohesive force	c	5 kPa
	Friction angle	$\Phi$	28°
	Bulk modulus	K	1.5 MPa
	Shear modulus	G	0.9 MPa
Tillage tool	Density	${\rho }_0$	2800 kg/m3
	Shear modulus	G	26 GPa
	Poisson’s ratio	$\upsilon$	0.33
	Melting temperature	$T_{melt}$	925 K
	Specific heat	$C_p$	875 J/(kg K)
	J–C plasticity model	A	324 MPa
		B	114 MPa
		N	0.002
		C	0.42
		m	1.34

, the representative properties of the soil and aluminum tool enable comprehensive studies of the tool’s performance and soil interactions under typical loading conditions relevant to agricultural applications. By incorporating a deformable tool and characterizing the progressive damage, this model provides enhanced capabilities beyond the traditional rigid tool assumptions used to study tool fracture failures.

3.2. Process of Tool Deformation and Fracture

This study utilizes numerical simulations to investigate the progressive development of deformation and the eventual ductile fracture failure of agricultural tillage tools under representative loading conditions. A virtual cultivation environment is established consisting of a simplified tillage tool cutting through a rectangular soil domain at a constant speed. As depicted in Figure 3b, the upper section of the tool is modeled as a non-deformable body that moves horizontally from right to left at a defined speed of 2 m/s. The lower cutting section of the tool is deformable and follows the motion of the rigid upper section. This setup aims to provide enhanced insights into the full progressive failure process, from initial elastic response through extensive plastic deformation to eventual fracture.

Initially, the front tip of the tool experiences minor elastic strains due to the cutting resistance concentrated at the leading end. As the tool advances further into the soil, the contact region between the tool and the soil gradually grows larger, resulting in the development of substantial plastic deformation along the lower portion of the tillage tool. The plastic zone expands to encompass the entire tool length, at which point the tool enters a plastic strengthening stage. Here, the deformation sharply escalates as the tool is subjected to increasing loads over larger contact areas. Ultimately, the excessive accumulated plastic strain exceeds the ductility limit of the tool material, leading to crack initiation and propagation. The virtual model enables close examination of how the ductile fracture process evolves from localized elastic strains to extensive plastic flows, culminating in catastrophic failure crack formation. By simulating the comprehensive interactive behavior between the deformable tool and soil medium, the model captures the continuous evolution of the tool’s deformation and damage response under representative loading conditions relevant to agricultural applications. The results illustrate the mechanisms leading to fracture failure from initial cutting resistance, offering an effective platform to study agricultural tool fracture and optimize tool designs.

Figure 4 provides a detailed visualization of the progressive development of deformation and eventual ductile fracture failure of the agricultural tillage tool throughout the simulated dynamic cutting process. Figure 4a shows the initial undeformed configuration of the tool as it first makes contact with the pristine soil medium prior to any cutting. The sequence in Figure 4b through to Figure 4h then demonstrates the continuous evolution of the tool’s deformation response over time as it interacts with the soil during the advancing cutting motion. These simulation results enable extensive analysis of the full progressive damage process leading to fracture of the tool. Specifically, Figure 4b reveals that minor localized elastic strain is first concentrated at the front tip of the tool at t = 2 ms after the initial soil contact. As cutting continues, Figure 4c through to Figure 4h show the tool deforming extensively due to the increasing soil resistance forces along the cutting interface. The accumulated soil build-up applies escalating loads along the surface of the tool. This leads to a transition from the initial localized elastic response at the tip to a substantial plastic flow extending across much of the tool length. The large plastic deformation culminates in ductile overload fracture, as shown in the final Figure 4h. By capturing the comprehensive interactive behavior, these detailed simulation visualizations provide unique insights into the complex transition from elastic deformation to extensive ductile damage. The model enables a thorough examination of the mechanisms by which increasing soil cutting forces first lead to reversible elastic strains, followed by the progressive development of severe plastic deformation that ultimately precipitates catastrophic fracture failure. This understanding of the complete failure sequence under representative loading facilitates the optimization of agricultural tillage tools to improve robustness and fracture resistance.

The deformation at the front tip of the tool is more obvious, and Figure 5 illustrates the deformation pattern here. In the early stage of tool deformation, the soil and tool contact area is small, and the cutting resistance is also small, so the front tip of the tool is not significantly deformed. As the contact area increases, the cutting resistance grows and the deformation at the front tip of the tool is similar to a linear change before fracture. The deformation pattern can provide a more accurate supplement as well as a reference for improving the robustness and fracture resistance of agricultural tillage tools.

To further investigate the tool’s resistance to deformation under loading, this study conducts comparative simulations utilizing tool models with varying mechanical properties. Specifically, the deformation response is analyzed for tools modeled with both aluminum alloy and steel materials exposed to identical simulated cutting conditions. The simulation results clearly demonstrate that the steel tool exhibits a substantially enhanced resistance to plastic deformation compared to the relatively softer aluminum tool, as shown by the deformation patterns in Figure 6. Under the same applied loads, the aluminum tool experiences more extensive plastic flow, while the steel tool is better able to withstand the forces with higher strength and hardness, resulting in reduced strains.

Quantitatively, the maximum deformation in the aluminum tool reached over 4.8 mm before failure, while the steel tool only deformed by approximately 1.9 mm plastically. This difference highlights the importance of material strength and hardness in avoiding excessive plastic damage and improving robustness against deformation and fracture failure. This comparative study provides an effective platform for evaluating tool materials and optimizing material selection to withstand anticipated cutting loads during agricultural tilling operations. Harder tool materials with higher yield strengths, such as wear-resistant steels, are better able to resist the soil cutting forces and exhibit enhanced fracture resistance compared to softer metals such as aluminum alloys.

The contrasting deformation patterns and strain levels further elucidate the mechanisms of fracture and damage progression in tools under cutting loads. While both materials accumulate strain over time, the limited ductility and earlier fracture of the aluminum alloy show it is more prone to damage accumulation and sudden unstable failure. In comparison, the steel tool undergoes lower progressive straining, allowing more uniform deformation distribution and avoidance of localized damage concentration. This highlights the substantial benefits of higher-strength steels in delaying the onset of fracture through greater damage tolerance. Overall, these fundamental insights into the relationships between tool material properties and deformation resistance will aid in the design and material selection of agricultural tillage tools. Selecting tool materials with adequate hardnesses, strengths, and fracture toughnesses is crucial to improving durability, damage resistance, and reliability under the harsh loading conditions inherent to agricultural tilling operations.

3.3. Distribution of Soil Displacement

The dynamic cutting motion of agricultural tillage tools induces complex extrusion and shear deformations in the soil, and shear deformations tend to occur at locations where the cumulative plastic strain is large. As shown in Figure 7, the formation of soil shear zones can be predicted by observing the changes in cumulative plastic strain.

When the soil is disturbed, the discrete soil particles are displaced in a highly non-uniform manner, depending on their spatial location. As depicted in the particle visualizations in Figure 8, the soil particles directly along the cutting edge experience the most intense compressive and shearing forces as the leading edge of the tool actively drives the soil flow. In contrast, the soil particles located below and along the length of the tool undergo lower overall displacements.

This is particularly evident after a fracture occurs, where the fractured segment of the tool substantially disrupts the underlying soil movement, causing decreased soil displacements in the middle and lower portions. This differential motion reflects the soil deformation field generated by the complex progressive interactions with the deforming and fracturing tool. Thoroughly analyzing these heterogeneous soil displacements provides the opportunity to relate localized soil deformation behaviors to the performance and fracture failures of the specific tool geometries and operating conditions. These fundamental insights can then be leveraged to optimize the tool designs and tillage parameters to achieve more uniform displacement distributions throughout the soil, potentially improving tillage effectiveness and efficiency. For instance, modifying tool geometries and cutting angles to reduce fractures may assist in producing more consistent soil displacements. By enabling detailed analysis of the complex coupling effects and spatial variations, this study facilitates the enhanced design of agricultural tillage tools for improved performance, fracture resistance, and tillage effectiveness.

3.4. Velocity Vector Analysis

Analyzing the velocity variations of the tool and the soil helps to verify the above results for the tool deformation and soil displacement distribution. Figure 9 shows the velocity vector plots of the tool and soil, where the color and length of the arrow line segments represent the velocity magnitude of each particle, and the direction of the arrow represents the velocity direction of each particle. The figure directly visualizes the velocity state of tool particles and soil particles at different moments.

The results show that the trend of tool deformation and soil displacement distribution are consistent with the velocity of the tool particles and soil particles. After the tool contacted the soil, the velocity of the lower portion of the tool particles changed, and the larger the velocity variable, the larger the tool deformation; meanwhile, the deformation trend was consistent with the velocity direction. The greater the velocity of the soil particles along the above cutting edge, the greater the soil particle displacement; the location of the soil accumulation formation was also consistent with the velocity direction.

3.5. Cutting Force Analysis

Cutting force is an important parameter that directly governs energy consumption, efficiency, and reliability during agricultural tillage operations. Accurately modeling the evolution of the resulting cutting forces on the tool over time is therefore essential for optimizing the performance and durability of tillage tools. In this study, the evolution of the cutting forces on the tool is simulated over time by modeling the complex interactive forces between the dynamic tool and discrete soil particles. Figure 10 presents the detailed simulated profiles of both the normal and tangential cutting force components developing along the tool surface under three different tool cutting speeds: 0.5 m/s, 1 m/s, and 2 m/s, respectively. The cutting force magnitudes and durations are shown to vary quantitatively based on the different tool speeds, reflecting differences in soil deformation rates and tool damage progression. Tracking the force development throughout the cutting process at these representative speeds enables in-depth studies of the relationships between tool geometry, operating parameters, and resulting loads on the tool. This facilitates the enhanced design of agricultural tillage tools for optimal energy efficiency and prolonged service lifetimes under the harsh soil-cutting environment parameters.

As shown through the cutting force profile comparisons in Figure 10, the overall qualitative trends and force evolution stages are relatively consistent across the different tool speeds analyzed. However, the force magnitudes and durations vary quantitatively based on the distinct cutting velocities. Examining the representative case with a cutting speed of 2 m/s provides detailed insights for correlating the progressive cutting force development with the tool’s deformation and damage processes. Initially, at 0.7 ms after the tool contacts the soil, the cutting force rapidly increases as the resisting soil applies increasing loads on the advancing tool. This sharp rising force trend is attributed to the growing soil cutting resistance as the tool–soil interface expands laterally and soil builds up along the tool surface. Around 2 ms, the onset of substantial plastic deformation in the tool structure causes a slight drop and fluctuations in the cutting forces. In the subsequent plastic strengthening stage, the cutting forces gradually increase again as additional loads are applied over larger contact areas. The forces peak at 7 ms as tool fracture initiates, followed by a declining trend as loads are relieved by the fracturing tool. Tracking this intimately coupled force–deformation history highlights the direct relationships between cutting forces and accumulating damage in the tool.

Combined with the trends of cutting forces during the tool deformation process and the deformation of the tool in the above conclusions, it is found that these results coincide with the four stages of metal deformation (elastic, yield, strengthening, and necking stages), thereby proving the correctness of the coupled soil–tool model proposed in this paper.

This understanding of these tight correlations facilitates the optimization of agricultural tillage tool designs for enhanced fracture resistance and durability. For instance, modifying geometries to sustain cutting loads while reducing damage accumulation will improve tool performance and longevity over prolonged operation cycles. Overall, these simulations capturing the interactive cutting physics provide an effective platform for studying design enhancements that can meet the demanding mechanical requirements of soil tilling equipment and operations.

4. Conclusions

This study utilizes the SPH method to simulate the progressive plastic deformation and ductile fracture failure of agricultural tillage tools. The coupled soil–tool model calculates the complex interactions and predicts the damage process:

(1)

The model combines an elastoplastic constitutive model and the Johnson–Cook damage criterion, incorporating Lennard-Jones repulsive forces for the soil–tool interface coupling. Numerical techniques including the artificial stress terms and Jaumann rate correction handle the challenges of cohesive soil, stress fluctuations, and large deformations. This enables simulation of the complete progressive failure sequence from initial loading to final fracture.

(2)

Simulations reveal the detailed post-fracture displacement fields in the soil, highlighting the non-uniform distributions. Velocity vector plots visualize and accurately reflect the motion of the tool and soil particles. The coupled force–deformation tool response provides insights into the relationships between cutting forces and accumulating tool damage. This fundamental understanding facilitates the optimization of tool structural design for enhanced fracture resistance.

(3)

Overall, the model accurately characterizes the complex soil–tool interactions involving substantial stresses and deformations. This technique could be extended to simulate progressive damage in other earthmoving equipment components, providing a versatile simulation platform.