Discrete Element Model Building and Optimization of Tomato Stalks at Harvest

Abstract

The mechanical properties of tomato stalk, relevant to the harvesting and crushing of tomato vines, significantly impact its harvesting quality and efficiency. Establishing a simulation model, which accurately mirrors these properties, is foundational for designing related mechanical components. The discrete element method models tomato stalk harvesting and is optimized through mechanical tests and simulations. A blend of Plackett–Burman, steepest ascent, and central composite design modeling identified three contact model parameters influencing the maximum stalk shear force. The optimal values of these three parameters were a normal stiffness of 1.04 × 10¹⁰ N m⁻³, tangential stiffness of 7.59 × 10⁹ N m⁻³, and bond radius of 1.06 mm. The relative error in the simulated versus measured shear force was <1%, affirming the model’s accuracy in characterizing cutting properties. These findings lay the theoretical groundwork for numerical simulations of tomato-stalk-related equipment.

1. Introduction

Tomatoes are among the three primary globally traded vegetables, together with onions and potatoes, and play a significant role in the global vegetable market. In 2020, tomato cultivation covered 5.055 million hectares worldwide, yielding 182.05 million tons. China is the world’s top tomato producer. In 2021, China cultivated tomatoes on 1.113 million hectares, producing 66.09 million tons, accounting for one-third of the global tomato production [1]. Tomato vines, the byproducts of tomato harvesting, serve as valuable biomass resources. After harvesting, these vines lie flat on the ground, characterized by their long, dense, and tangled forms. Tomato stalks possess intricate mechanical properties and high strength. When crushed by pulverization, they offer substantial resistance to the rotating blades. Optimizing the crushing efficiency and uniformity of tomato stalks is essential to enhance resource utilization [2]. Therefore, studying the mechanical properties of tomato stalks aids in the design of straw harvesters and provides crucial theoretical data.

In recent years, the discrete element method (DEM) has been widely used to calibrate crop simulation parameters and reveal internal mechanical properties, which has been driven by advancements in computer technology. Wei et al. (2023) used EDEM software to create a DEM model for the double-layer bonding of rapeseed shoot stalks. They refined it through mechanical testing and demonstrated its suitability for stalks and similar crops [3]. Tong et al. (2023) calibrated the parameters for a DEM model of corn stover using EDEM, highlighting the role of simulation software in analyzing crop mechanical properties [4]. Vahid and Chen (2018) developed a DEM-based model for the tensile properties of industrial hemp fibers using indicators such as tensile strength and Young’s modulus for calibration [5]. Kovács (2019) used the Timoshenko beam-bond model to create a DEM model for corn stalks and investigated their compression, shear, and bending properties [6]. These studies do not include a methodology for modeling discrete elements of tomato stalks, nor do they provide a complete set of discrete element parameters as a reference; however, they do provide insights, which can be used for conducting this research.

In this study, greenhouse tomato stalk fractures during shear were simulated using the DEM with harvested stalks. Modeled as particles with homogeneous material properties, tomato stalks were simulated using EDEM software. The mechanical parameters were calibrated through measured and simulated tests to explore the mechanical properties of stalks and enhance model reliability. This study aims to improve the numerical simulation accuracy of tomato stalk force dynamics during harvesting, providing a theoretical basis for simulating similar materials and machinery designs for stalk collection and crushing.

2. Materials and Methods

On 12 April 2023, a random sample of greenhouse cherry tomato vines was collected from a vegetable production base on Guli Street, Jiangning District, Nanjing City, Jiangsu Province, China. Healthy, fresh plants without pests or diseases were carefully selected, and efforts were made to minimize stalk damage. The morphological structure of a single tomato vine plant is depicted in Figure 1, comprising leaves, lateral branches, the main stalk, and roots. The main stalk was subcircular in the cross-section. Subsequently, the samples were harvested, with the leaves and lateral branches removed, leaving the middle stalk as the primary focus for testing.

The tests used a WDW-20-type microcomputer-controlled electronic universal testing machine from Jinan Chuan Bai Instrument Co., Ltd, Jinan, China. Tomato stalks underwent shear and bending tests. The specimens were 80 mm long and obtained from the middle section of the stalks to ensure uniformity and integrity. Over 20 sets of stalks were prepared for the shear and bending tests. In the shear tests, ten 80 mm tomato stalks were selected and numbered. The test points were located at the midpoints. A uniform downward load was applied at a rate of 25 mm min⁻¹. Carbon steel inserts with a thickness of 15 mm and an edge angle of 30°, as shown in Figure 2a, were used. A three-point bending method was applied to measure the elastic and shear moduli of the stalks in the bending tests. The outer diameter and bast thickness were measured at the center of each stalk. A continuous vertical load of 1 mm min⁻¹ was applied during the bending tests, as illustrated in Figure 2b.

The shear force–displacement values were derived directly from the shear tests. In the bending test, after obtaining the compressive force–displacement values, the modulus of elasticity and shear modulus were calculated. The cross-section of the stem was approximated as circular, and the moment of inertia concerning the neutral axis was calculated using Equation (1). The elastic modulus of the stalk was determined using Equation (2). The Poisson ratio of the tomato stalk was assumed to be 0.2, based on previous studies [7,8,9,10]. The shear modulus of the stalk was calculated using Equation (3). The test data obtained from the shear and bending tests are presented in

Table 1. Basic physical parameters and mechanical properties of tomato stalks.

Parameter	Average Value	Standard Deviation
Stem length of whole plant (mm)	2090	3.2
Stalk diameter (mm)	11.1	0.63
Stalk bast thickness (mm)	3.64	0.56
Core diameter (mm)	4.92	0.49
Stalk moisture content (%)	75.9	4.2
Stalk density (kg m−3)	769.96	23
Maximum shear force (N)	220.57	3.52
Elastic modulus of stalk (MPa)	50.64	2.42
Shear modulus of stalk (MPa)	21.09	2.81

.

$$I=\frac{\pi \left[d_1^4-{\left(d_1-2\omega \right)}^4\right]}{64}$$

(1)

$$E_{\omega }=\frac{F{L}_1^3}{48SI}$$

(2)

$$G=\frac{E}{2\left(\mu +1\right)}$$

(3)

where I is the moment of inertia of the stalk cross-section around the center axis; E_ω is the elastic modulus; G is the shear modulus; d₁ is the stalk outer diameter; ω is the stalk bast thickness; F is the loading force; L₁ is the distance between the two supports; S is the bending deflection at the midpoint of the stalk; and μ is the Poisson ratio of the stalk.

The DEM, which was introduced by Cundall and Strack, is a numerical simulation technique [11]. It involves a system composed of numerous particles, with time divided into discrete steps. During each time step, a four-step process is executed comprising neighborhood search, contact determination, contact force analysis, and numerical integration. This approach effectively simulates the system changes over time [12]. In the EDEM 2022 version of the simulation software, the Hertz–Mindlin (no slip) contact model, standard rolling friction model, and Bonding V2 bonding model were used to simulate tomato stalk movement, damage, and fracture. The interaction forces between the particles were determined using the Hertz–Mindlin (no slip) model, which serves as the default model in EDEM, as depicted in Figure 3. The mathematical expressions for the Hertz–Mindlin (no slip) model with respect to the interaction forces are as follows:

$$F_n=\frac{4}{3}{E}^{*}\sqrt{R^{*}}{\delta }_n^{\frac{3}{2}}$$

(4)

$$F_n^d=-2\sqrt{\frac{5}{6}}\beta \sqrt{S_n{m}^{*}}{v}_n^{\overrightarrow{rel}}$$

(5)

$$F_t=-8G^{*}\sqrt{R^{*}{\delta }_n}{\delta }_t$$

(6)

$$F_t^d=-2\sqrt{\frac{5}{6}}\beta \sqrt{S_t{m}^{*}}{v}_t^{\overrightarrow{rel}}$$

(7)

where F_n is the normal force; $F_n^d$ is the normal damping force; F_t is the tangential force; $F_t^d$ is the tangential damping force; E* is the equivalent Young’s modulus; R* is the equivalent radius; δ_n is the normal overlap; δ_t is the tangential overlap; m* is the equivalent mass; $v_n^{rel}$ is the normal component of the relative velocity; S_n is the normal stiffness; G^* is the equivalent shear modulus; and $v_t^{\overrightarrow{rel}}$ is the relative tangential velocity.

In the Bonding V2 model, a virtual beam cell is formed between two spherical particles (Figure 4), denoted as the bond. The bond transfers the forces and moments between the particles. When the stress exceeds a predetermined strength, the bond fails and cannot be reintroduced [11]. The contact radius should be set to 20–30% above the actual physical radius of the particles [13]. The bond breaks when the stress exceeds the maximum value.

$${\delta }_{max}<\frac{-F_n}{A}+\frac{2M_t}{J}{R}_b$$

(8)

$${\tau }_{max}<\frac{-F_t}{A}+\frac{M_n}{J}{R}_b$$

(9)

where δ_max is the maximum normal stress; τ_max is the maximum tangential stress; R_b is the bond radius; M_t is the normal shear stress; and M_n is the tangential shear stress.

The tomato stalk simulation model was constructed using basic spherical particles in the EDEM simulation software. An appropriate particle-contact radius was set when creating the spherical particles. Bonds were established between these particles, imparting a degree of flexibility to the discrete stalk meta-model. To ensure accuracy, the model considers the actual mechanization and essential physical properties of the tomato stalk, with references considered. The determinations of Poisson’s ratio, density, and shear modulus for both tomato stalks and steel in this test were based on prior research [10,14,15,16,17,18]. The intrinsic parameters are detailed in

Table 2. Intrinsic and surface contact parameters of tomato stalk.

Symbol	Value
Poisson’s ratio of stalk	0.2
Density of stalk (kg m−3)	769.960
Shear modulus of stalk (MPa)	21.098
Poisson’s ratio of steel	0.3
Density of steel (kg m−3)	7800
Shear modulus of steel (MPa)	7.00 × 104

.

The physical dimensions of the particles were 0.75 mm for the radius and 0.9 mm for the contact radius, as determined from actual measurements of the tomato stalk. This decision was aimed at enhancing simulation efficiency and reducing computational time. To maintain a small model error, we minimized the number of particles while ensuring an accurate representation of the desired physical structure.

This was achieved by approximating the tomato stalk to a cylinder. The cross-sectional plane of the stalk was set as the xy-plane, which served as the particle generation plane within the right-handed Cartesian coordinate system. The z-axis represented the axial direction of the stalk, with the origin O (0, 0, 0) denoting the center coordinates of the cross-sectional circle. The position of each particle was determined by the sphere center’s x-, y-, and z-coordinate values. Cross-sectional particles were distributed in two-layered circles in non-nodal regions and in three-layered circles at stalk nodes. All constructed particle cross-sections were meticulously aligned along the axial direction without any gaps, resulting in the creation of the desired physical structure.

Customized functions and the Excel 2023 software autofill feature were used for the computation of all 3D coordinate data, adhering to the above-described method. Prior to the incorporation of the meta-particle function, we set parameters such as Poisson’s ratio, shear modulus, and density for the individual spherical particles in the EDEM software. We generated a complete simulation model of a tomato stalk by importing point data from a coordinate data table. Figure 5a,b illustrate the simulation model of the stem shear specimen. Subsequently, the contact parameters for the bond key were added based on this model, as depicted in Figure 5c,d, as the bond composition in the stem simulation model.

In the stem mechanics simulation, the corresponding stem model had the same length as the actual specimen. In the simulation calculations, the auxiliary fixed components, which were irrelevant to the stem mechanics tests, were removed. The three-dimensional software SolidWorks 2022 was used to establish the corresponding test fixture and shear tool model, which were then imported into the EDEM software. Simulation accuracy and efficiency were carefully balanced to ensure the continuity of the simulation. A fixed time step of 1.57 × 10⁻⁶ s was set for the simulation calculations. Based on the measured test method, the tool was loaded onto the stem simulation model at a loading rate of 25 mm min⁻¹ with a preload force of 0.2 N. The discrete element geometry used for the stem mechanics simulations is shown in Figure 6.

This test used the maximum shear as the parameter calibration index. The simulation parameters were analyzed and calibrated using a combination of physical measurements and simulations. Ranges of values for the contact parameters were selected to perform the stem simulation shear tests. The simulated shear force was derived from the vertical force applied to the shear knife in the mechanical simulation model, with the maximum shear force as the target value for the test. Simulation approximations and comparative analyses were performed. The Plackett–Burman test, the steepest ascent method (SAM), and the central composite design (CCD) were sequentially used to optimally adjust the parameters of stem simulation [19]. The test evaluation index was the minimum relative error between the simulated and actual maximum shear forces.

3. Results

The tomato stalk simulation model used in this study required the specification of several contact and bonding parameters during the construction process in EDEM. These parameters included the coefficient of restitution between stalk and stalk, coefficient of restitution between steel and stalk, coefficient of static friction, coefficient of rolling, and bond parameters. Initially, the value ranges of these parameters were derived from simulated parameters of similar materials, such as vines and cereal straws, as referenced from previous studies [9,20,21,22,23]. Subsequently, several simulated shear pretests were conducted, and the parameter ranges were appropriately expanded. The final values of the contact and bonding parameters were set within the ranges shown in

Table 3. Contact parameters of tomato stalk DEM model.

Parameter	Symbol	Parameter Levels
		−1	1
Coefficient of restitution between stalk and stalk	A	0.24	0.36
Coefficient of static friction between stalk and stalk	B	0.56	0.84
Coefficient of rolling friction between stalk and stalk	C	0.08	0.12
Coefficient of restitution between steel and stalk	D	0.24	0.36
Coefficient of static friction between steel and stalk	E	0.40	0.60
Coefficient of rolling friction between steel and stalk	F	0.08	0.12
Normal stiffness per unit area (N m−3)	G	9.60 × 109	1.44 × 1010
Shear stiffness per unit area (N m−3)	H	8.00 × 109	1.20 × 1010
Normal strength (Pa)	I	1.92 × 107	2.88 × 107
Shear strength (Pa)	J	1.60 × 107	2.40 × 107
Bonded disk	K	0.80	1.20

.

Minitab 21 software facilitated the design of the Plackett–Burman screening test, taking maximum shear as the response variable. The differences between the high and low levels of the test factors were compared to the total difference. The contact and bond parameters underwent screening to determine their effect on the shear force. In the Plackett–Burman design, the high levels of each factor do not typically exceed twice the low level [24]. The design scheme and simulation results of the Plackett–Burman design are shown in Table 3, with parameters A–K representing the coded values of each parameter.

Drawing on the Plackett–Burman design, the results are displayed in

Table 4. Scheme and results of the Plackett–Burman design.

No.	A	B	C	D	E	F	G	H	I	J	K	Y (N)
1	1	−1	1	1	−1	−1	−1	−1	1	−1	1	241.91
2	1	1	−1	1	1	−1	−1	−1	−1	1	−1	87.63
3	−1	1	1	−1	1	1	−1	−1	−1	−1	1	222.93
4	−1	−1	1	1	−1	1	1	−1	−1	−1	−1	174.26
5	1	−1	−1	1	1	−1	1	1	−1	−1	−1	131.96
6	1	1	−1	−1	1	1	−1	1	1	−1	−1	75.12
7	1	1	1	−1	−1	1	1	−1	1	1	−1	138.38
8	1	1	1	1	−1	−1	1	1	−1	1	1	300.66
9	−1	1	1	1	1	−1	−1	1	1	−1	1	187.02
10	1	−1	1	1	1	1	−1	−1	1	1	−1	84.85
11	−1	1	−1	1	1	1	1	−1	−1	1	1	400.29
12	1	−1	1	−1	1	1	1	1	−1	−1	1	231.48
13	−1	1	−1	1	−1	1	1	1	1	−1	−1	150.51
14	−1	−1	1	−1	1	−1	1	1	1	1	−1	131.17
15	−1	−1	−1	1	−1	1	−1	1	1	1	1	183.65
16	−1	−1	−1	−1	1	−1	1	−1	1	1	1	393.21
17	1	−1	−1	−1	−1	1	−1	1	−1	1	1	160.00
18	1	1	−1	−1	−1	−1	1	−1	1	−1	1	279.74
19	−1	1	1	−1	−1	−1	−1	1	−1	1	−1	94.14
20	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	−1	104.02

. A regression model significance analysis for the maximum shear was performed, and the results are shown in

Table 5. ANOVA for the Plackett–Burman design impact rates.

Source	df	Adj SS	F-Value	p-Value
Model	11	163,671	10.65	0.001 **
A	1	4789	3.43	0.101
B	1	499	0.36	0.567
C	1	1270	0.91	0.368
D	1	633	0.45	0.520
E	1	701	0.50	0.499
F	1	845	0.60	0.459
G	1	39,640	28.38	0.001 **
H	1	11,593	8.30	0.020 *
I	1	87	0.06	0.809
J	1	1532	1.10	0.326
K	1	102,083	73.08	0.000 **
Pure error	8	11,175
Cor total	19	174,847

Note: ** indicates a highly significant effect (p < 0.01), and * indicates a significant effect (p < 0.05).

. The regression model for maximum shear, Y, is expressed in Equation (10):

Y = 188.65 − 15.47A + 5.00B − 7.97C + 5.63D + 5.92E − 6.50F + 44.52G − 24.08H − 2.09I + 8.75J + 71.44K

(10)

The analysis of variance (ANOVA), shown in Table 5, and the pareto plot (Figure 7) underscore the model’s significance with p < 0.01 and R² = 0.94, indicating a significant main effect. This suggests that the regression equation modeled is consistent with the experiment. The factors from A to K have a variable influence on the maximum shear response. Factors K, G, and H are highly significant, while others are less so. As shown by the regression Equation (1), factors K and G exert positive effects, and factor H exerts a negative effect. Therefore, calibration and optimization were exclusively performed for the three significant bonding parameters, K, G, and H, in subsequent SAM and CCD procedures. The remaining non-significant factors were set at their intermediate levels, namely A = 0.3, B = 0.7, C = 0.1, D = 0.3, E = 0.5, F = 0.1, I = 2.4 × 10⁷ Pa, and J = 2 × 10⁷ Pa.

Significant parameters were screened using the Plackett–Burman design. SAM was used to quickly approach the optimal vicinity. The relative error (ξ) between the maximum shear obtained through SAM (F_a) and the average maximum shear from the real test (F_b = 220.57 N) was used as the criterion for the optimal solution in the final CCD test.

$$\xi =\frac{\left|F_a-F_b\right|}{F_b}\times 100\%$$

(11)

Table 6. Scheme and results of the steepest ascent method.

No.	G	H	K	Fa (N)	ξ (%)
1	9.60 × 109	8.00 × 109	0.80	90.74	58.86
2	9.98 × 109	7.82 × 109	0.88	126.84	42.49
3	1.04 × 1010	7.63 × 109	0.96	184.39	16.41
4	1.07 × 1010	7.45 × 109	1.04	229.03	3.83
5	1.11 × 1010	7.27 × 109	1.12	271.88	23.26
6	1.15 × 1010	7.08 × 109	1.20	299.43	35.75

shows the SAM design scheme and results. The results showed that the smallest relative error of 3.83% occurred at test level No. 4. To accurately derive the tomato stalk contact parameters, CCD was conducted using the SAM parameter from Group 4 as the central value, with experimental data from Groups 3 and 5 providing the range of level values. The maximum shear was used as the response variable to seek the optimal solution.

Following the SAM results, the tomato stalk contact parameters were further explored through CCD to derive the maximum shear force and its relative error. The effects of the CCD maximum shear on the normal stiffness, tangential stiffness, and bond radius among the tomato stalks were examined.

Table 7. Factor coding for CCD simulation tests.

Levels	Factors
	G	H	K
−1.68	1.0364 × 1010	7.265 × 109	0.96
−1	1.05 × 1010	7.34 × 109	0.99
0	1.07 × 1010	7.45 × 109	1.04
1	1.10 × 1010	7.56 × 109	1.09
1.68	1.1036 × 1010	7.635 × 109	1.12

shows the coding for the CCD simulation test factors. The CCD design scheme and its results are shown in

Table 8. Programs and results of CCD.

No.	G	H	K	Fmax (N)	ξ (%)
1	1	1	−1	201.35	8.71
2	−1	−1	−1	202.61	8.14
3	0	1.68	0	220.84	0.12
4	1	−1	1	255.82	15.98
5	0	0	1.68	256.77	16.41
6	0	−1.68	0	234.82	6.46
7	0	0	0	229.78	4.18
8	0	0	0	229.78	4.18
9	0	0	0	229.78	4.18
10	1	1	1	252.01	14.25
11	−1	1	−1	199.09	9.74
12	−1	−1	1	241.21	9.36
13	0	0	−1.68	191.69	13.09
14	0	0	0	229.78	4.18
15	−1.68	0	0	216.09	2.03
16	0	0	0	229.78	4.18
17	1.68	0	0	216.09	2.03
18	−1	1	1	237.58	7.71
19	1	−1	−1	218.16	1.09
20	0	0	0	229.78	4.18

.

The experimental results underwent binary regression analysis, which is detailed in Table 8. A second-order regression model for the maximum shear force was developed, correlating with the coded values of three factors:

F = 229.66 + 3.43G − 3.75H + 20.12K − 1.68G·H + 1.40G·K + 1.61H·K − 4.06G² + 0.09H² − 1.18K²

(12)

Table 9. ANOVA of the CCD test.

Source	Sum of Squares	df	Mean Square	F-Value	p-Value
Model	6196.15	9	688.46	36.99	<0.0001 **
G	160.81	1	160.81	8.64	0.0148 *
H	192.51	1	192.51	10.34	0.0092 **
K	5532.26	1	5532.26	297.22	<0.0001 **
GH	22.70	1	22.7	1.22	0.2953
GK	15.78	1	15.78	0.8476	0.3789
HK	20.78	1	20.78	1.12	0.3155
G2	237.78	1	237.78	12.77	0.0051 **
H2	0.1168	1	0.1168	0.0063	0.9384
K2	20.20	1	20.2	1.09	0.3221
Residual	186.14	10	18.61
Lack of fit	186.14	5	37.23
Pure error	0.0000	5	0
Cor total	6382.29	19

Note: ** indicates a highly significant effect (p < 0.01), and * indicates a significant effect (p < 0.05).

shows the CCD results, which were analyzed by ANOVA. The results show a second-order regression model with p < 0.0001, a coefficient of determination R² = 0.9708, and a corrected coefficient of determination R² adj = 0.9446, both of which were close to 1. This shows a highly significant relationship between the dependent variable (maximum shear) and all independent variables. The order of each factor’s influence on the maximum shear was C > B > A. Only the A² interaction term had a significant effect on the maximum shear force. The effects of the remaining interaction terms were statistically insignificant.

Using design-expert 13 software, we obtained the response surfaces for the interactions between each factor and the maximum shear force from their respective regression equations, as depicted in Figure 8. Specifically, Figure 8a illustrates the response surfaces of tangential stiffness and bond radius at a constant normal stiffness. The figure shows that the maximum shear force incrementally increased with the bond radius when the tangential stiffness remained constant. However, at a constant bond radius, the impact of an increase in tangential stiffness on the maximum shear force was less pronounced.

Figure 8b illustrates the response surfaces for normal stiffness and bond radius at a steady tangential stiffness. The maximum shear force was observed to rise with an increase in bond radius at a consistent normal stiffness. Conversely, at a fixed bond radius, an increase in normal stiffness led to a gradual increase in the maximum shear force.

Figure 8c presents the response surfaces for normal and tangential stiffnesses at a fixed bond radius. Here, with constant normal stiffness, the maximum shear force decreased as the tangential stiffness increased. However, holding the tangential stiffness constant, the maximum shear force initially increased and then decreased with rising normal stiffness, although this trend was not highly pronounced.

Based on the CCD results and the quadratic regression equation, the objective was to minimize the relative error between the simulated and the actual average maximum shear forces. An optimal solution was sought for each factor, with the objective function and constraints set as follows:

$$\left\{\begin{array}{c}min \xi \left(G,H,K\right)\\ s.t.\left\{\begin{array}{c}-1.68\le G\le 1.68\\ -1.68\le H\le 1.68\\ -1.68\le K\le 1.68\end{array}\right.\end{array}\right.$$

(13)

Finally, one set of optimized solutions with relative errors < 1% was obtained. The coded values are (G = −1.665, H = 1.235, and K = 0.443), corresponding to the simulation parameters of normal stiffness at 1.0367 × 10¹⁰, tangential stiffness at 7.5859 × 10⁹, and bond radius at 1.0611. These were used as the optimal contact parameters. Simulated shear tests conducted with these optimal contact parameter combinations yielded a curve delineating the relationship between shear force F and radial strain of the stalk, which was then analyzed and compared with real test results. The test results showed that at all stages of shearing (beginning, middle, and end), phenomena such as bond compression by shear, neat cut-off, and cross-sectional rupture concurred with observed real-life phenomenon, thereby providing a better characterization of tomato stalk shear damage morphology. The tomato stalk shearing process with the bond breaking process shown in Figure 9 was obtained. Figure 10 shows a comparative analysis between the measured mechanical test results and simulation test results for the radial strain–shear force:

$$\epsilon =\frac{\mathsf{\Delta }b}{b}$$

(14)

where Δb is the shear displacement, and b is the stem diameter.

A comparative analysis of the F–ɛ plot (Figure 10) revealed that the overall trends of the simulated and measured curves align consistently. This relationship was linear during the rising phase of the curve, and the curve experienced a rapid decline after reaching the peak shear value. Oscillations in shear values within this range were evident in both the measured and simulated data. This suggests that at that point, both the measured and simulated stalks have cross-sectional breaks. The average measured maximum shear force was 220.57 N, while the maximum shear force of the simulated curve was 220.84 N, resulting in a relative error of 0.12%. These results show that the optimally combined contact parameters can accurately reflect the shear mechanical properties of stalks, confirming the simulation method as feasible and highly reliable.

4. Discussion

In this study, we developed a discrete elemental model of the tomato stalk using EDEM software, aiming for high particle accuracy [10,25]. The flexibility characteristics, eigen parameters, and contact parameters of our model are similar to those found in previous studies [11,12,19,23]. Tomato stalks in previous studies have been viewed as solid stalks, but tomato stalks differ from corn and rapeseed stalks; they lack a full medulla and have a hollow, approximately circular cross-section. In our modeling process, we considered the moment of inertia of the cross-section around the central axis to ensure that our simulation parameters corresponded closely with the actual mechanical properties of tomato stalks.

In our calibration study, which focused on the DEM-based flexible contact parameters of the stalk, we concentrated on understanding the variations in normal stiffness, tangential stiffness, and bond radius parameters. Our study was limited to the middle stalks of tomato plants at harvest time, using the shear loading method to investigate the contact parameters. While this approach has its limitations, it serves as a foundational step for our ongoing research. Future experiments will extend to a comparative analysis of different components of the tomato vine plant, including the root, upper, lower, and leaf materials, with a focus on varying water content. Additionally, we plan to enhance DEM modeling by separately characterizing and calibrating the bast and medulla under different loading conditions, including tension, compression, bending, and shear, to provide a more comprehensive representation of the entire tomato vine plant.

5. Conclusions

• We initiated this study by establishing the fundamental physical parameters of the tomato stalk, creating a discrete meta-model, which closely approximated real stalks, and constructing a geometric model to simulate stalk mechanics, mirroring an actual shear-test setup.
• Our evaluation criterion was the relative error between the simulated maximum shear force and the observed average maximum shear force in the physical tests. We used simulation approximation and comparative analysis techniques to calibrate the simulation parameters for the tomato stalks. The Plackett–Burman design underscored the significance of three factors—normal stiffness, tangential stiffness, and bond radius—for the maximum shear force. For precise parameter calibration, we conducted the SAM within specific ranges for normal stiffness, tangential stiffness, and bond radius: 1.03 × 10¹⁰–1.11 × 10¹⁰ N m⁻³, 7.27 × 10⁹–7.63 × 10⁹ N m⁻³, and 0.96–1.12 mm, respectively. Using CCD and ANOVA, we obtained a second-order regression model for these three factors. After setting the objective functions and constraints, we analyzed the optimal solutions for each factor. The process yielded solutions with relative errors below 1% and determined the following parameter values: normal stiffness at 1.04 × 10¹⁰ N m⁻³, tangential stiffness at 7.59 × 10⁹ N m⁻³, and bond radius at 1.06 mm.
• In this study, we used the EDEM simulation software to simulate the discrete element contact model for tomato stalks. The simulation results closely align with the measured curves, reflecting the trend and maximum values accurately. Furthermore, we successfully replicated the rupture characteristics of stem ligaments and the flexibility observed in actual experiments.

This research significantly enhances the accuracy of numerical simulations depicting the mechanical stress dynamics of tomato stalks during mechanical harvesting. The outcomes of our investigation provide valuable theoretical insights for the future development of numerical simulations and the design of tomato stalk collection and crushing equipment.