Kaysse Ibrahim ^1,2*; Mayssa Shash ¹; Jalal Abboud ²

1, Department of equipment and machinery engineering, Faculty of Technical Engineering, University of Tartous, Tartous, Syria

2, General Commission for Scientific Agriculture Research (GCSAR), Hejaz Station, Damascus, Syria

 

Manuscript link
http://dx.doi.org/10.30493/DAS.2024.485466

Abstract

This study introduces a comprehensive dynamic model for abrasive wear of the threshing tooth. The ensuing contact-impact issue was articulated and resolved as a linear complementarity problem, integrated within Moreau’s time-stepping method. Local contact parameters, including normal contact forces, tangential stresses, and wear distribution for each contact surface cell, were estimated utilizing Archard’s wear model. The wear progression was simulated in MATLAB and executed in a loop, with the worn tooth profile being updated in the dynamic model after every 106 load cycles. The impact of the threshing cylinder’s rotational velocity and the friction coefficient was examined, and the operational lifespan of the rasp-bar was assessed. Predictions of tooth wear were analyzed and contrasted with empirical facts. The wear prediction model, as detailed further in this work, can ascertain the optimal timing for re-profiling or replacement well in advance of the critical wear depth specified in the maintenance standards.

Keywords: Threshing unit, Wear, Dynamics, Multibody systems

Introduction

Wheat, a crucial strategic crop globally, requires meticulous focus on its growth, harvesting, and post-harvest procedures including threshing, packaging, and storage. The threshing process is essential for maintaining grain quality and reducing losses. Despite improvements in threshing machine designs, the deterioration of essential components, especially rasp-bars, continues to pose a substantial challenge, adversely impacting machine performance, efficiency, and the quality of the harvested grain.

The design of agricultural gear is fundamentally dependent on comprehending the physical and mechanical characteristics of grains. Previous research has examined the impact of factors including friction coefficients, moisture content, and surface rigidity on the performance of threshing machines [1–5]. Nevertheless, these studies frequently emphasize the optimization of machine settings while inadequately considering the wear phenomena, which gradually affects operating efficiency and threshing quality over time.

The effectiveness of threshing is influenced by various elements, including the torque exerted on the threshing cylinder, friction and impact forces, and the mechanical interaction between the crop material and machine components. Excessive friction and impact forces result in grain damage and energy losses while also hastening the wear of rasp-bars. Previous studies have demonstrated that grain damage is significantly associated with pre-impact velocity and the stiffness of contact surfaces [6–8]. Nonetheless, the impact of prolonged usage on these forces and overall machine efficacy has garnered little scrutiny.

This problem transcends technological issues and presents considerable economic difficulties. The deterioration of threshing components escalates operational expenses due to increased energy consumption and the necessity for frequent maintenance or replacement. Additionally, wear-induced inefficiencies can result in increased grain losses and diminished grain quality, especially for grains designated for seed utilization. Moreover, increased wear imposes further stress on machines, amplifying fuel and energy use while diminishing total productivity. Consequently, mitigating wear is crucial for reducing operational expenses, preserving grain quality, and enhancing economic results [9][10].

Researchers have created multiple computational models to examine the mechanical behavior of grains subjected to threshing pressures. Finite element methods have been employed to examine stress distributions and mechanical properties under diverse situations [11–15]. The impact of dynamic vibrations on threshing teeth was previously invistigated [16], whilst Pfister and Eberhard analyzed frictional contact issues in elastic and stiff systems [17]. Glocker and Pfeiffer formulated models for multi-impact dynamics [18], whereas Lizhang et al. emphasized the significant impact velocity influencing grain damage [19]. Nevertheless, these studies frequently regard threshing as discrete impact events, overlooking the ongoing friction interactions between crop material and machine components. This omission neglects to account for the wear dynamics of rasp-bars and their influence on tooth-concave clearance, a crucial factor for effective threshing [20][21].

The degradation of rasp-bars not only enhances tooth-concave clearance but also modifies the geometry of threshing teeth, hence diminishing the overall efficacy of the threshing unit. Calibration is rendered useless when the tooth height is less than the acceptable threshold (≥6 mm). This highlights the necessity for extensive research on the wear phenomena and its impact on machine performance and grain quality.

The present study addresses this gap by introducing a dynamic model for the abrasive wear of threshing teeth. Utilizing the Linear Complementarity Problem (LCP) framework and Moreau’s time-stepping method, this research simulates the contact-impact interactions during threshing and estimates wear distribution using Archard’s wear model. The findings provide insights into wear dynamics and their implications for machine performance, offering a foundation for improving threshing machine designs and operational strategies, while also mitigating economic losses associated with inefficiencies and grain damage

Material and Methods

Frictional impact kinematics

Threshing transpires due to the friction between the wheat ears and the rasp-bars as the crop mass traverses the gap between the threshing cylinder and the concave, where σ denotes the concave winding angle and DC represents the concave diameter, as illustrated in (Fig. 1 A) [22]. The clearance value varies from di= 16 mm at the crop entry point to de= 8 mm at the exit position [23] (Fig. 1 B).

In order to describe the location and direction of the rasp-bar tooth and the crop mass with respect to the coordinate system, an OXY inertial coordinate system was built to represents the generalized coordinates of the system and the thresher frame oxy at the point o (Fig. 2).

To simplify analyses, it was assumed that at a certain velocity there was an impact between the plant mass and the rasp-bar with a short-duration slip motion. Considering ψ the rotation angle of the threshing cylinder, and θ =(Ψ+Π/2), the coordinates of the potential contact point p on the surface of the rasp-bar tooth in the oxy coordinate system can be expressed in h_p (x_p,y_p). r_p is the location of the potential contact point on the rasp-bar on the oxy frame. The coordinates in the OXY frame of the potential contact point p, and the rotational transformation matrix (W_P)from the OXY to the oxy frame, can be calculated as follows:

where g_Tk is the tangential gap and g_Nk is the normal gap in the base coordinates, and thereafter:

In which r_Q is the location of the potential contact point on the crop mass on the O’X’Y’ frame (Fig. 3) W_Q is the rotational transformation matrix from the OXY coordinate system to O’X’Y’ frame, (k) is the potential contact pair between the rasp-bar tooth and the crop mass, l the distance between the reference coordinate center o and the generalized coordinate center O, n^T_Pk and τ^T_Pk are the normal and tangential vectors of the contact pair k, respectively.

By substituting Eq. 1, 2, and 4 into Eq. 3, the tangential gap  and the normal  in the base coordinates can be obtained:

By differentiating Eq. 5 and assuming that the vegetation mass is static and motionless, Ẇ_Q=0, ṙ_q=0, the relative velocities in the tangent and normal directions of the contact pair can be obtained:

For potential collision, and to describe the various collision states of contact pairs, four index metrics are defined as follows [24]:

I_G = {1,2,……,n_G}

I_S = {k ∈ I_G∣g_nk = 0,k =1,2,….,n_S}

I_C =  {k ∈ I_S∣ġ_nk = 0,k =1,2,….,n_C}

I_H = {k ∈ I_C∣ġ_tk = 0,k =1,2,….,n_H}

The set I_G represents all the points n_G in a single collision. The set I_S indicates that the normal gap collision point n_S is zero, and the tangential relative velocity is arbitrary value. The set I_C indicates that the normal constraint n_C is in a constant and continuous contact state including all sliding and stick states, in which the normal relative velocity is zero. The group I_H indicates that the contact points n_H are in the state of stick/slip switch, n is the number of elements that can change during each time step.

The crop inside the threshing unit takes the form of a tape with a thickness equal to the clearance between threshing cylinder and the concave (δ) (Fig. 4) [25] with the same way as mentioned before, W_Q.r_Q can be obtained by,

where (x_q,y_q) is the position of O’ with respect to the generalized coordinates OXY, and the angle φ is the rotational transformation angle, which its value in this case is zero where always (O’X’ // OX).

The frictional path in the threshing unit (the path of the potential contact point on the rasp-bar surface when the threshing phase in active) is a circular path with center O and radius l_i. (where, i= 1,2…,26), which can be written with respect to the generalized coordinates as follows,

where x_i and y_i are the generalized coordinates of the potential impact point on the rasp-bar tooth when it meets the thresher concave.

The profile of the threshing rasp-bar tooth is a quadric equation with three constants (Fig. 5) is estimated as follows:

In this case of study, the values of the constants are (α=-5.539, β=-0.1014, γ=0.0188), where the tooth height is equal to 1.88 cm and the tooth width is equal to 5 cm.

Formalisms for the rasp-bar undergoing large overall rotation

Considering han arbitrary position of the center of mass of the rasp bar in the frame OXY, J the moment of inertia of the rasp-bar, m the mass of the rasp bar, l_m the distance between the center of mass and the center of the generalized coordinates O, (θ ̇=ω)  is the angular velocity of the rasp-bar with respect to the OXY coordinates. The system kinetic energy T is defined as follows:

The system gravitational potential energy can be obtained by,

where h_g is the height of the position of the center of gravity from the lowest position in the dynamic system, g is the gravity acceleration (9.81 m/s²), the behavior of the dynamic system of the separation phase is described by the Lagrange method:

Where L is the Lagrangian dynamic system, by deriving Eq. 16 for the generalized coordinates θ:

By taking the second derivative, generalized momentum can be obtained by:

The system behavior of the separation state is generally described by the Lagrange method:

where R is the generalized force. The gyroscopic and external forces are given by:

Where h: the total generalized forces acting on the rasp-bar during the separation phase. : the operative torque.

Formalisms for dynamics of the frictional impact of the multibody system

The time-stepping methods provide a discrete numerical scheme suitable for non-smooth systems simulation [26-28]. These methods are widely used due to their simplicity to implement and their robustness. The time-stepping schemes are based on a time-discretization of the system dynamics. The whole set of discretized equations and constraints is used to compute the next state of the motion.

The dynamic equations of motion of a multi-body systems (MBS) with normal and tangential contact forces during an impact can be written at the acceleration level as [29]:

where  is the positive definite and symmetric mass matrix,  and  represent the generalized normal and tangential force directions, respectively. The normal and tangential contact forces have magnitudes λ_Ti and λ_Ni for each contact point i. u= addresses the system generalized velocities (equal to ω) and ụ is the vector that contains the system accelerations.

where t_M is the midpoint time instant of the compact time interval [t_A, t_E] and q_M is the midpoint system’s position state. The midpoint time instant can be evaluated as:

The equations of motion expressed at the velocity level and in a small finite time interval Δt, of which A and E denote the initial and end points [27], can be written a:

Where M_M=M(q_M,t_M), h_M=h(q_M,u_A,t_M), w_NM=h(q_M,t_M), w_TM=h(q_M,t_M) are the matrix notation that were used above, μ is the coefficient of friction, and H is the set of active contacts at moment t,u_E, and u_A are vectors of the generalized velocities at t_E and t_A, respectively.

This set of algebraic inclusions can be solved with a linear complementarity problem (LCP) formulation by augmented Lagrangian approach (ALA). Then the velocity u_E, at the end of time-step t_E=t_A+Δt, is subsequently calculated using Eq. 23. Finally, the positions at the end of the time step are calculated:

Considering unilateral contacts, Moreau’s midpoint algorithm calculates the contact distances g_Ni of all unilateral contacts at the midpoint q_M to evaluate whether they are active (g_Ni≤0) or not (g_Ni>0). Only active unilateral contacts can be modelled by inclusion (24). Unilateral contacts that are non-active, thus open, are disregarded because it is assumed that their contact force contribution is equal to 0.

At present, there is no consensus on how to determine the initial conditions for collisions [21]. In this study, it was assumed that the normal and tangential restitution coefficients are equal to zero, ε_n=0 and ε_t=0, respectively; This is to ensure continuous contact between the rasp-bar tooth and the mass of the flowing crop, which completes the conversion from non-collision to collision. The impulse-momentum method is derived from the classical rigid-body collision theory, which is a method for solving collisions based on the generalized momentum equation and the restitution coefficient equation.

The normal and tangential impulse arrays of the collision point are  And . Let the normal victor n_k and the tangent τ_k of the contact pair k be consistent with the multi-body, respectively. For all impacts in the system k∊I_S, n = diag (n₁,n₂,……, n_ns) and τ = diag(τ₁,τ₂,……,τ_ns), where (s) represents the number of contact points in a single impact.

The jump time is infinitesimal, therefore, the tangential and radial collision impulse, respectively P_T and P_N of the impact point can be considered as a constant force in the context of the impact. The following two equations give the normal and tangential dynamic impulsive indicators according to [29]:

The contact-impact problem of the multibody system is solved to determine the impact and contact forces using the linear complementarity problem (LCP), which follows to form the appropriate work algorithm along with the wear size prediction equations according to Archard and is embedded in the Moreau’s time-stepping method [27],  LCP formulation is presented to solve the contact-impact problem for multibody systems with unilateral frictional constraints [27].

In order to set up the LCP, the following matrix notation, can be described as follows,

The contact-impact problem of non-smooth systems can be summarized as follows:

The relative normal and tangential contact velocities  and  can be calculated using Eq. 31 and Eq. 32, respectively, since the velocities u_A are known at the left endpoint of the time interval. By introducing Eq. 29 and 30 into 33 and 34:

Where ξ_N and ξ_T are two defined parameters to account for the possibilities of the contact does not participate in the impact, that is, the value of the normal contact impulse is zero, although the contact is closed. This situation happens normally for multiple contact scenarios. Therefore, for this case, it is allowed that the post-impact relative velocity to be higher than the value prescribed by Newton’s impact law, with the intent to express that the contact is superfluous and could be removed without changing the contact-impact process [30]. w ̃_N and w ̃_T are the Jacobian terms representing the rheonomic constraints and they are here equal to zero, since the constraints are not rheonomic.

The inclusions for the contact-impact force laws can be formulated as complementarity conditions. Thus, the unilateral primitive of Eq. 24 results in:

In turn, the relay function (35) must be decomposed into two Uprs to achieve the desired complementarity conditions:

in which the step height is [-μP_N,+ μP_N]. Additionally, to abbreviate the complementarity conditions of Eq. 39 the impulsive friction saturations P_R and P_L are defined as [31]:

It should be noted the special arrangement of P_L and ξ_L in Eq. 42, and they must be placed in this manner, which has deep roots in optimization theory. Otherwise, the LCP formulation cannot be set up without additional matrix inversion processes [31]. Since variables ξ_L, P_T and u_E are not included in (42), they must be eliminated. Thus, combining Eq. 28 and Eq. 40 yields:

The elimination of variable P_T can be done through the combination of Eq. 40 and Eq. 41, and it can be written as:

Since the inversion of mass matrix M is always possible, then Eq. 42 can be solved for u_E:

Now, Eq. 31 and Eq. 32 are used to express W_NM^Tu_A and W_TM^Tu_A in terms of γ_NA and γ_TA :

Introducing Eq. 45, Eq. 46 and Eq. 47 into Eq. 36 and 43 yields:

Thus, Eq. 48, Eq. 49 and Eq. 44 can be written in a matrix form as:

Eq. 50 together with the complementarity conditions (Eq. 42) form the LCP for the contact-impact analysis of multibody systems with frictional unilateral constraints. The dimension of this LCP is 3i, where i represents the number of active contacts. The LCP (Eq. 50) is solved at each integration time step. Then, the velocities u_E and positions q_E for the subsequent time steps are obtained from Eq. 45 and Eq. 26, respectively.

It should be highlighted that for Δt=0, the LCP reduces to the pure impact equations of motion and can be used, for example, for initialization of the velocities. For = =0, the LCP describes impact free motion at the velocity level, containing still the cases of persisting contact and stiction as well as transitions to sliding or separation and can be transformed to the acceleration level by the methods presented by [31].

Rasp-bar tooth wear estimation

Archard’s wear model is considered a suitable solution in abrasive wear applications to depict the surface wear phenomenon since it depends on a simple linear relationship between the volume of material loss, contact pressure and slip distance. More precisely, this model is based on a phenomenological approximation of the abrasive wear mechanism that considers the applied load and the slip distance between the contacting bodies. Archard’s model was chosen in this study due to its compatibility with the dynamic model of the non-smooth multi-body system. This presented model considers wearing constant over the time. Therefore, the Archard’s wear equation can be expressed as follows [32]:

where V is wearing volume (m³), K is wearing coefficient depending on the contact conditions and contact pairs (m²/N), F_N the applied load (N), s is the slip distance (m). In this study, Archard’s wear equation was introduced into the MATLAB environment and combined with the algorithm for the dynamics of the non-smooth multibody system. For a representative elementary volume of material, Eq. 51 can be written as:

where dV is elemental wear volume, ds is the elemental change of the slip distance (Fig. 6). The elemental wear depth (dW) and contact pressure (p) can be clearly inferred by introducing the elemental contact area d𝛺 as follows:

The wear depth for each contact point can be written as follows:

The contact area between the threshing rasp-bar and the mass of the flowing grains increases with increasing the normal force to reach a maximum value of 3×10^-8 m² [31].

The simulation provides data at each time step [t, t+Δt], therefore, the requested variables were written in predefined frequency for computational efficiency, so the simulation was divided into several frames, The n step corresponds to one frame and the step (n+1) is related to the following frame:

where,

The sliding distance during each time step s_piⁿ can be calculated as follows:

where γ_TAi is the tangential velocity, and ∆t is the simulation time step it expresses the sliding distance during each time step, which can be inferred by determining the slip velocity during each step for one step time.

The wear depth is calculated at each time step, and this value is added periodically according to the Euler integration method. After calculating the amount of worn material, the tooth profile must to be updated and smoothened for further dynamic analysis after numerous load cycles (equal to the number of revolutions of the threshing cylinder per 1 h), thus, updating the values of α, β, γ. The wear estimation along with the updating tooth profiles is carried out again into the algorithm, where the y_p value is renewed until the final number of simulation cycles is reached according to the following equation:

where y_i is the original tooth height at the contact point p for all cells of the tooth profile i. Since the number 200 indicates the number of working hours, given that a single cycle is equal to a single working hour, to finally estimate the wear depth, the value of the wear coefficient (K) was 5×10^-17 m²/N, which was practically calculated by measuring the wear rates of the rasp-bar by the Archard’s wear model after 13 working hours by calculating the weight loss for the rasp-bar before and after the test.

To compare simulated results with experimental observation, a threshing unit was used is a tangential feeding and axial flow thresher that was specially designed and manufactured for this study by The General Commission for Scientific Agricultural Research (GCSAR), Syria (Fig. 7). The feed rate was previously calibrated to be 5.2 kg/s, the output torque of the thresher and the speed were detected simultaneously using a torque-speed sensor. the detected torque on the applied feed rate and rotation speed was 146.49±6.2 Nm and 1145.8±13 rpm, respectively, the threshing was performed on winter wheat with an average height of 0.87 m and the moisture content of the stem and grain were 13% and 12%, respectively. The rasp-bar was removed and cleaned from rust and dust, washed with water and dried. Then, a digital camera (Optilia w30x) was used to acquire digital images of tooth profile before and after the test.

Computational strategy for solving kinematic equations and estimating wear depth

The dynamic model of the threshing process was carried out based on the presented algorithm and created in MATLAB environment. Then, the model was embedded in the Moreau’s time-stepping method during the total simulation time t_f with a small-time step (10^-6 s). The angular velocities of the threshing cylinder were 1145.8, 954.86, 859.37, 740, 763.8, and 668.4 rpm (i.e. 120, 110, 100, 90, 80, 70 rad/s). The effective torque was 146.49, 110.91, 54.37, and 20.87 Nm for each angular velocity. The rasp-bar inertia was 0.0005186444 kg.m², the weight of the rasp-bar was 0.346 kg, and the coefficient of friction was 0.14, 0.18, 0.25, 0.35, 0.3, 0.4.

The value of t_F was 0.02 s, which is a sufficient time to complete a single load cycle. The initial angle of movement was ψ = 15^o, considering that the surface roughness of the thresher tooth is equal to zero, and there are no penetrations between the contact surfaces, the flow rate of the crop mass was assumed to be constant.

Since Moreau’s time-stepping method with LCP formulation includes a good set of mathematical formulations, it is convenient to summaries the main steps with a convenient algorithm (Fig. 8), that was created under the multi-body systems formulation. To solve the kinematic equations for solid multi-body systems with unilateral frictional constraints [34], the wear depth/distribution was calculated by the Archard’s wear model. The algorithm can be summarized by the following steps:

1. Enter all the data such as the coordinates of the surface points under study, which are 26 points/cells, the constants  of the rasp-bar tooth section equation α, β, γ, the initial angle ψ, the contact surface area d𝛺, the gravitational acceleration g, and the number of load cycles per minute that corresponds to the number of the threshing cylinder revolutions per minute.
2. Start the analysis by defining the initial conditions, start moment of the Moreau’s step interval t_a, final moment of one load cycle t_f, step time Δt, initial position q_A and initial velocity u_a.
3. Define the geometrical, inertia and material functions g_Ni, gyroscopic and external forces vector h_M, the normal and tangent restitution coefficient ε_Ni, ε_Ti, the friction coefficient μ_i, the dynamic normal and tangential impulsive indicators w_Ni, w_Ti, the Jacobian terms representing the rheonomic constraints w ̃_Ni and w ̃_Ti, and the value of the wear coefficient K. The initial number of load cycles is equal to 1, and the initial wear depth is equal to 0.
4. Compute the frictional impact kinetics of non-smooth solids with unilateral constraints.
5. According to the Moreau’s midpoint rule, compute the midpoint time moment t_m, end time of the interval t_E, evaluate the position at the midpoint moment q_M, prepare the mass matrix MM for the midpoint and vectors of gyroscopic and external forces h_M (Eq. 20), compute the midpoint states variables of all potential contact points H_M;       

6. Check for impact-contact, If there is no contact-impact (open contact) the velocity is calculated at the end of time, u_E, by using the Eq. (23).
7. Otherwise (at least one closed contact) apply the time-stepping with (LCP) algorithm to obtain the normal P_N and tangent P_T impulse forces needed to calculate the final velocity u_E of the contact-impact state and the volume loss, this procedure is performed on each cell of the twenty-six cells on surface of the tooth.
8. For each:
   

9. Setup LCP with the standard form of y_LCP=A.x_LCP+b_LCP usingEq. 50.
10. Solve LCP using the appropriate algorithm: (x_LCP,y_LCP)=LCP(A,b).
11. split the LCP solution according to:

12. Calculate wear depth w_pi and the tooth height  using Eq. 57-60.
13. Velocity evaluation at the end of the integration time step using Eq. 45.
14. Calculate the locations at the end of the integration time step using Eq. 26.
15. Time step increment, if the current time is less than the final simulation time, update the position and velocity parameters, then go to step (6) to complete the new time step process; Otherwise the simulation is stopped: t_A=t_A+∆t.
16. Increase the number of load cycles, then calculate the wear depth at the end of the load cycle and update the coordinates of the studied points.
17. Update the values of α, β, γ and smoothening the equation of the rasp-bar tooth section, and the system states’ variables q_A=q_E and u_A=u_E, and then move to step (2).
18. This loop is repeated until the number of maximum load cycles is reached (cycleF=200).

Results and Discussion

The influence of rotational speed of the threshing unit

The wear rate of the rasp-bar was analyzed after 4×10⁶ load cycles at varying rotational speeds and applied torques. The results indicate a linear increase in wear rate with rotational velocity up to 90 rad/s for 110.39 Nm and 100 rad/s for 145.39 Nm (Fig. 9). Beyond these points, wear rates rise significantly, with higher torques amplifying the effect. At lower speeds, the differences in wear are minimal, but they become more pronounced as speed increases. This trend is driven by the combined effects of tangential and normal sliding velocities, which elevate impulse forces and wear rates.

At high torques and rotational speeds exceeding 90 rad/s, the likelihood of slip conditions increases, resulting in higher contact stresses and accelerated wear. This behavior aligns with findings from prior studies. One of these studies showed higher energy consumption, wear, and grain breakage at increased rotor speeds [35]. Another work demonstrated that enhanced threshing performance at higher speeds comes at the cost of increased wear [36]. Additionally, it was noted that higher rotor speeds improved separation efficiency but also increased wear and energy consumption [37], mirroring the trends observed in the current study.

Wear distribution along the tooth length (Fig. 10) varied, with maximum wear observed at reference point 6 and minimal wear at point 26. This uneven wear pattern is due to localized frictional stresses and interaction with the plant mass. These findings align with other studies, which highlighted the influence of forces distribution along a threshing drum on accelerated wear in regions of high material contact and stress​ [38].

The largest wear rates appeared at the highest rotational speeds where they reached 7620, 8280, 9020, 9940, 13400, and 17070 μm, at rotational velocities of 70, 80, 90, 100, 110, and 120 rad/s, respectively at 146.39 Nm, and reached 1130, 1410, 1640, 2050, 2500, and 2670 μm, respectively at 20.87 Nm. Hence, the shape and curvature of solid surfaces with unilateral frictional constraints have an important effect on the effective momentum forces and wear rates [39].

The worn profiles with the number of applied load cycles at each rotational speed at a torque 146.39 Nm are given in Fig. 11. From the simulated data, it appears that the wear mechanism is the same at all velocities, however, the wear depth becomes greater with the velocity increases as explained previously, and the tooth profile begins to change at higher rotational speeds (110-120 rad/s), when wear is stronger. The negative number of the tooth height means that the tooth has been completely worn out and the wear starts to appear at the base of the rasp-bar, which has a maximum thickness of 8 mm. From these results it is possible to predict the effective life of the rasp-bar at a specific rotational speed.

The influence of the coefficient of friction

The influence of the coefficient of friction (μ) on wear rates was analyzed across various torque levels, revealing a clear correlation between increasing μ and the acceleration of wear. At lower torques (20.87–54.37 Nm), wear occurred more uniformly along the tooth surface, with wear depth ranging between 492–718 μm at μ = 0.14 and 2569–2850 μm at μ = 0.4. As torque increased, wear patterns became more localized, with significant wear near areas of higher curvature, such as reference point 6 (Fig. 12). This is consistent with findings that higher contact stresses in curved regions amplify wear intensity. At higher torques (e.g., 146.39 Nm), the wear depth increased dramatically with increased μ values, confirming that curvature amplifies stress concentration. This observation aligns with other studies on agricultural machinery, where higher contact stresses in curved components enhance abrasive interactions, accelerating localized wear​ [40].

The influence of the friction coefficient (μ = 0.4) on the wear value appears clearly with increasing load cycles, as its highest values appear at the reference point 6 (Fig. 13). However, at the end of the simulation, the concentration of this wear decreases at point 6, which in turn confirms that the greater the curvature of the friction surfaces, the higher the value of the stresses as well, and thus, the greater the wear values.

When the coefficient of friction increases above 0.18, the intensity of the wear increases (Fig. 14). This means that the lower the crop moisture, the higher the wear values, because the coefficient of friction is directly related to the moisture of the crop [4].

By converting the load cycle into working hours, and taking the rotational velocity into account, the life of the rasp-bar can be estimated. For example, when μ = 0.14, the efficiency of the rasp-bar ends after 3.5 million load cycles (64 working hours) and the rasp-bar becomes completely ineffective after 3.8 million load cycles (69.5 working hours) (Fig. 13). When μ = 0.18, the tooth efficiency ends after 3.6 million load cycle. When μ = 0.25, the efficiency of the rasp-bar ends after 1.5 million cycles, i.e, after 27 working hours and must be replaced after 40.5 working hours. When μ = 0.3, the rasp-bar efficiency ends after 31 working hours, and after only 14 working hours when μ = 0.35, and 20 working hours when μ = 0.4. After these limits, the rasp-bar wears out faster, and it can be noticed that the rasp-bar wears out completely at the end of the simulation for μ = 0.35 and μ = 0.4.

Comparison with the experimental results

The simulated tooth profiles were evaluated after 2×10⁶, 3×10⁶, and 4×10⁶ load cycles (at a rotational speed of 70 rad/s) and after 1×10⁶, 2×10⁶, and 3×10⁶ load cycles (at a rotational speed of 120 rad/s). Then, simulated results were compared with experimental data (Table 1). The results show that the mean relative error (ep​) increases with the number of load cycles. The greatest discrepancies between simulated and experimental tooth heights were observed at the front and middle of the tooth for 70 rad/s and 120 rad/s, respectively (Fig. 15). These discrepancies can be attributed to the changing influence of tooth geometry on wear depth due to variations in unit vector directions at different rotational speeds. The maximum mean relative errors were 8.928% and 8.604% for 70 rad/s and 120 rad/s, respectively.

The observed differences between simulated and experimental results can be explained by several factors:

1. Variations in the coefficient of friction, which may fluctuate as wear progresses during testing.
2. Experimental inconsistencies in rotational speed, applied torque, and feed rate.
3. Measurement errors introduced during the conversion of tooth profile images into numerical data.
4. The potential presence of fretting wear, which was not considered in the simulation based on Archard’s wear model, as this model assumes linear wear evolution.

despite these discrepancies, the simulation results provide a reliable prediction of the abrasive wear phenomena in rasp-bar teeth. The experimentally derived worn tooth profile exhibits a cubic curve, with larger wear depths in specific regions, whereas the simulation used quadratic approximations (Fig. 15).

It is important to note that external environmental conditions, such as temperature and humidity, also influence wear behavior. For example, elevated temperatures can soften materials, increasing wear rates, while high humidity may act as a lubricant, reducing wear through decreased friction. Although these factors were not explicitly analyzed in this study, they warrant consideration in practical applications and future research.

Conclusions

This study presented a comprehensive model to predict wear behavior in the rasp-bar tooth during the threshing process, solving the contact-impact problem using a linear complementarity framework embedded in Moreau’s time-stepping method and utilizing Archard’s wear model for wear distribution analysis. The model successfully estimated wear depth across the contact surface, offering valuable insights into optimal re-profiling or replacement schedules before reaching critical wear depth thresholds specified in maintenance guidelines.

The results highlight key findings:

1. Rotational Speed Influence: Wear rates increase linearly with rotational speed and torque up to specific thresholds (90 rad/s for 110.39 Nm and 100 rad/s for 145.39 Nm). Beyond these speeds, wear accelerates dramatically due to higher tangential and normal sliding contact velocities.
2. Coefficient of Friction: The coefficient of friction exerts a dominant effect on wear compared to rotational speed. Higher friction coefficients, typically associated with reduced crop moisture, significantly amplify wear depth, especially at high torque levels.
3. Wear Distribution: Uneven wear distribution across the tooth length was observed, with maximum wear at the reference point nearest the leading edge. This aligns with localized frictional stress concentrations and the geometry of the rasp-bar.

The model effectively predicted rasp-bar life under different operating conditions. For example, at low friction coefficients (e.g., μ = 0.14), the rasp-bar remains functional for 64 working hours. In contrast, at higher friction coefficients (e.g., μ = 0.35), functional life drops to 31 working hours.

Comparative Analysis with Experimental Results: The simulated results closely matched experimental data, with mean relative errors between 2.01% and 8.93%, depending on the rotational speed and load cycles. Discrepancies were attributed to factors such as fluctuating friction coefficients during testing, potential measurement errors, and external environmental influences like temperature and humidity. Despite these limitations, the simulation reliably predicts abrasive wear patterns and life expectancy of rasp-bars under diverse conditions.

Practical Applications: The methodology enables:

• Optimizing rasp-bar design by comparing different geometries and materials under varying operational parameters.
• Extending machine availability by predicting maintenance needs, thus minimizing downtime and operational costs.
• Enhancing threshing efficiency by identifying optimal operating conditions that balance wear reduction and dynamic performance.

The algorithm’s versatility allows it to be adapted for analyzing different rasp-bar profiles and materials, aiding in the development of durable and efficient threshing systems. These findings pave the way for advancements in agricultural machinery design, improving both performance and sustainability in crop processing operations.

Acknowledgments        

The authors are grateful for the support of The General Commission for Scientific Agricultural Research – Syrian Arab Republic.

References

1. Kalkan F, Kara M. Handling, frictional and technological properties of wheat as affected by moisture content and cultivar. Powder Technol. 2011;213(1–3):116–22. DOI
2. Kashaninejad M, Ahmadi M, Daraei A, Chabra D. Handling and frictional characteristics of soybean as a function of moisture content and variety. Powder Technol. 2008;188(1):1–8. DOI
3. Shafaei SM, Kamgar S. A comprehensive investigation on static and dynamic friction coefficients of wheat grain with the adoption of statistical analysis. J. Adv. Res. 2017;8(4):351–61. DOI
4. Al-Mahasneh MA, Rababah TM. Effect of moisture content on some physical properties of green wheat. J. Food Eng. 2007;79(4):1467–73. DOI
5. Shafaei S, Heydari A, Masoumi A, Sadeghi M. Determining and modeling of static friction coefficient of some agricultural seeds. Jordan J. Agr. Sci. 2015;11(4).
6. Reyns P, Missotten B, Ramon H, De Baerdemaeker J. A review of combine sensors for precision farming. Precis. Agric. 2002;3:169-82. DOI
7. Marvin R. Paulsen. Fracture Resistance of Soybeans to Compressive Loading. Trans. ASAE. 1978;21(6):1210–6. DOI
8. Wang XR, Li YM, Xu LZ. Relationship between thresher velocities and rice grain broken rate. Trans. Chin. Soc. Agric. Eng. 2007;23(8):16-9.
9. Wang Y, Li D, Nie C, Gong P, Yang J, Hu Z, Li B, Ma M. Research Progress on the Wear Resistance of Key Components in Agricultural Machinery. Materials. 2023;16(24):7646. DOI
10. Yang S, Li W. The Impact of Socialized Agricultural Machinery Services on Land Productivity: Evidence from China. Agriculture. 2022;12(12):2072. DOI
11. Mohapatra G, Sahay SS. Wear and Tribology in Agricultural Machinery. In Friction, Lubrication, and Wear Technology. ASM International. 2017;984–1002. DOI
12. Li XP, Gao LX, Ma FL. Analysis of finite element method on mechanical properties of corn seed. Trans. CSAM. 2007;38(10):64-7.
13. Miu PI, Kutzbach HD. Modeling and simulation of grain threshing and separation in threshing units—Part I. Comput. Electron. Agric. 2008;60(1):96–104. DOI
14. Miu PI, Kutzbach HD. Modeling and simulation of grain threshing and separation in axial threshing units: Part II. Application to tangential feeding. Comput. Electron. Agric. 2008;60(1):105-9. DOI
15. Lizhang X, Yaoming L, Zheng M, Zhan Z, Chenghong W. Theoretical analysis and finite element simulation of a rice kernel obliquely impacted by a threshing tooth. Biosyst. Eng. 2013;114(2):146–56. DOI
16. Glocker Ch, Pfeiffer F. Dynamical systems with unilateral contacts. Nonlinear Dyn. 1992;3(4):245–59. DOI
17. .Fouda T, Bahnas O, Soltan M, Ghoname M. A study on threshing device wearing behavior for rice harvesting combine. Econ. Eng. Agric. Rural Dev. 2021;21(1):305-12.
18. Han I, Gilmore BJ. Multi-body impact motion with friction-analysis, simulation, and experimental validation. J. Mech. Des. 1993;115(3):412–22.
19. Pfister J, Eberhard P. Frictional contact of flexible and rigid bodies. Granular Matter. 2002;4(1):25–36. DOI
20. Glocker Ch, Pfeiffer F. Multiple impacts with friction in rigid multibody systems. Nonlinear Dyn. 1995;7(4):471–97. DOI
21. Qian Z, Jin C, Zhang D. Multiple frictional impact dynamics of threshing process between flexible tooth and grain kernel. Comput. Electron. Agric. 2017;141:276–85. DOI
22. Ivan G, Vladut V, Ganea-Christu I. Improving threshing system feeding of conventional combine harvesters. In International Symposium on Agricultural Engineering. 2015;43.
23. Sewell AJ. Some effects of concave to drum clearance and concave design on small grain threshing drum performance. J. Agric. Eng. Res. 1990;46:207–17. DOI
24. Leine RI, Glocker C. A set-valued force law for spatial Coulomb–Contensou friction. Eur. J. Mech. A/Solids. 2003;22(2):193–216. DOI
25. Steponavičius D, Kemzūraitė A, Kiniulis V, Zokaitis K. Analysis of corn ears movement within the threshing crescent using high-speed recording method. J. Meas. Eng. 2017;5(4):240-8.
26. Moreau JJ. Unilateral Contact and Dry Friction in Finite Freedom Dynamics. Nonsmooth Mechanics and Applications. Springer, Vienna. 1988;1–82. DOI
27. Glocker C, Studer C. Formulation and Preparation for Numerical Evaluation of Linear Complementarity Systems in Dynamics. Multibody Syst. Dyn. 2005;13(4):447–63. DOI
28. Pfeiffer F, Foerg M, Ulbrich H. Numerical aspects of non-smooth multibody dynamics. Comput. Methods Appl. Mech. Eng. 2006;195(50–51):6891–908. DOI
29. Pfeiffer F, Glocker C. Multibody dynamics with unilateral contacts. Springer Science & Business Media. 2000.
30. Glocker C. On frictionless impact models in rigid-body systems. Philos. Trans. R. Soc. Lond. A Math. Phys. Eng. Sci. 2001;359(1789):2385–404. DOI
31. Glocker C. Set-valued force laws: dynamics of non-smooth systems. Springer Science & Business Media. 2013. DOI
32. Hambli R. Blanking tool wear modeling using the finite element method. Int. J. Mach. Tools Manuf. 2001;41(12):1815–29. DOI
33. Horabik J, Molenda M. Force and contact area of wheat grain in friction. J. Agric. Eng. Res. 1988;41(1):33–42. DOI
34. Flores P, Leine R, Glocker C. Modeling and analysis of planar rigid multibody systems with translational clearance joints based on the non-smooth dynamics approach. Multibody Syst. Dyn. 2009;23(2):165–90. DOI
35. .Kiniulis E, Patarovas V,  Pranckeviciene I. Effect of the rotor speed on energy consumption and wear of threshing equipment. Agric. Eng. Int. CIGR J. 2014;16(2), 1-10.
36. Olaye AR, Moreira J, Hounhouigan J, Amponsah SK. Effect of threshing drum speed and crop weight on paddy grain quality in axial-flow thresher (ASI). J. Multidisciplin. Eng. Sci. Technol. 2016;3:3716-21.
37. Que K, Tang Z, Wang T, Su Z, Ding Z. Effects of Unbalanced Incentives on Threshing Drum Stability during Rice Threshing. Agriculture. 2024;14(5):777. DOI
38. Fu J, Xie G, Ji C, Wang W, Zhou Y, Zhang G, Zha X, Abdeen MA. Study on the distribution pattern of threshed mixture by drum-shape bar-tooth longitudinal axial flow threshing and separating device. Agriculture. 2021;11(8):756. DOI
39. Dotta ALB, Costa CA, Farias MCM. Tribological behavior of alumina obtained by low-pressure injection molding using factorial design. Tribol. Int. 2017;114:208–20. DOI
40. Lian C, Zhao W, Wang S, Zhang P, Dong X, Li Y. Friction and Wear Properties of Wheat Straw Powder Third Body on Steel-Steel Friction Pair. Agriculture. 2022;12(6):882. DOI

Cite this article:

Ibrahim, K., Shash, M., Abboud, J. Dynamic simulation of abrasive wear in grain threshing units. DYSONA – Applied Science, 2025;6(1): 200-222. doi: 10.30493/das.2024.485466