Establishment of Whole-Rice-Plant Model and Calibration of Characteristic Parameters Based on Segmented Hollow Stalks

Abstract

To address the limitations of the current discrete element model of rice plants in terms of accurately reflecting structural differences and threshing characteristics, this study proposes a whole-rice-plant modeling method based on segmented hollow stalks and establishes a whole-rice-plant model that accurately represents the bending and threshing characteristics of the actual rice plant. Initially, based on the characteristics of the rice plant, the rice stalk was segmented into three sections of hollow stalks with distinct structures—namely, the primary stalk, the secondary stalk, and the tertiary stalk—ensuring that the model closely resembles actual rice plants. Secondly, the mechanical and contact parameters for each structure of the rice plant were measured and calibrated through mechanical and contact tests. Subsequently, utilizing the Hertz–Mindlin contact model, a multi-dimensional element particle arrangement method was employed to establish a discrete element model of the entire rice plant. The bending characteristics of the stalk and the threshing characteristics of the rice were calibrated using three-point bending tests and impact threshing tests. The results indicated calibration errors in the bending resistance force of 4.46%, 3.95%, and 2.51% for the three-section stalk model, and the calibration error for the rice model’s threshing rate was 1.86%, which can accurately simulate the bending characteristics of the stalk and the threshing characteristics of the rice plant. Finally, the contact characteristics of the model were validated through a stack angle verification test, which revealed that the relative error of the stacking angle did not exceed 7.52%, confirming the accuracy of the contact characteristics of the rice plant model. The findings of this study provide foundational models and a theoretical basis for the simulation of and analytical applications related to rice threshing and cleaning.

1. Introduction

Rice is the third largest food crop globally and serves as a primary food source in China [1]. Currently, most regions in China have achieved near-complete mechanization of the rice harvesting process [2]. The interaction between rice plants and harvesting machinery directly impacts the quality of rice grains [3,4,5]. However, the closed structure of the threshing device makes it challenging to directly observe the movement of rice within the threshing chamber [6]. Additionally, the seasonal nature of rice cultivation imposes both temporal and financial constraints on the optimization of mechanical devices, significantly hindering research into the mechanisms of the harvesting process and the development of harvesting machinery [7]. Therefore, studying the biomechanical properties of rice and modeling the entire plant is crucial for overcoming the constraints related to the spatial and temporal dimensions of the harvesting process.

With the advancement of the discrete element method, its application in the field of agricultural engineering has become increasingly prevalent. Numerous scholars have investigated the modeling and application of various crops [8,9,10,11]. Unlike other plants, rice plants consist of stalks, spikes, leaves, and roots. The structure of the plant is complex, and the mechanical properties of each component vary significantly [12]. In terms of rice modeling research, Lensert and Leblicq et al. interconnected several rigid stalk models using spring-damping virtual keys to simulate the bending behavior of the stalk [13,14,15]. However, this modeling approach exhibited noticeable discontinuities during the simulation process, preventing the achievement of continuous and smooth bending deformation. Consequently, the researchers proposed increasing the number of bonding nodes to better approximate the bending deformation of the stalk model to that of real rice. Mao et al. proposed a rice stalk model that connects elastic hollow cylindrical keys between each spherical particle to simulate the hollow structure characteristics of rice stalks [16]. Shi et al. carried out mechanical tests on the stalk and established a stalk model characterized by anisotropic failure properties [17]. Schramm et al. proposed a method for calculating the local damping coefficient and bond Young’s modulus of a flexible stalk model, which enhanced the accuracy of model calibration [18]. Su and Liu et al. developed a hollow flexible stalk model suitable for threshing based on the Hertz–Mindlin model, which more accurately reflected the true characteristics of the rice stalk compared to the rigid stalk model [19,20,21]. Wang and Tang et al. investigated the effect of moisture content on the characteristics of the rice model and established discrete element models for grains and stalks under varying moisture contents [22,23]. Tang et al. utilized three-dimensional laser scanning to extract the outline of rice grains, subsequently creating a discrete element model of the grains through particle filling [24]. Similarly, other researchers have developed stalk models for banana, hemp, citrus, and cotton [25,26,27,28], as well as grain models for camellia oleifera, corn, and millet [29,30,31]. However, the aforementioned modeling methods only addressed specific parts of the plant structure and did not establish an whole-plant model. As a result, these models lacked overall coherence and were unsuitable for the simulation analyses of continuous deformation processes.

In the terms of whole-plant modeling research on rice, Chen et al. developed a simple yet flexible whole-plant model composed of rice stalks and rice grain clusters [32]. Liu et al. established a rice model using the method of particle polymerization, which can be used to study the continuous threshing process [33]. Xu et al. examined the physical characteristics of the entire rice plant and established an upright flexible rice plant model [34]. Wang et al. developed a model of the entire rice plant utilizing hollow cylindrical bonding structures and characterized the crushing behavior of flexible rice plants [35]. However, the aforementioned modeling methods differed from existing rice plant models and actual rice plants. The particle arrangement methods employed in these rice plant models were insufficient to simulate and analyze the separation of grains from the stalk. Therefore, the structural modeling and biological feature simulation of entire rice plants remained areas that required further exploration and refinement.

The purpose of this study is to establish an whole-rice-plant model that accurately reflects the bending and threshing characteristics of an actual rice plant, thereby enhancing the credibility of the simulation. Initially, this study measures the density, moisture content, contact characteristics, and mechanical properties of each structure of the rice plant, based on the characteristics of the rice plant. Subsequently, utilizing the Hertz–Mindlin contact model, a multi-dimensional element particle arrangement method is employed to establish a discrete element model of the entire rice plant. The bending and threshing characteristics of the model were calibrated using a three-point bending test and an impact threshing test. Finally, the contact characteristics of the model were validated using stack angle verification tests.

2. Materials and Methods

2.1. Test Material

This study utilizes the rice plant variety Shuyounuo 81, cultivated in Shijing Village (109.86 E, 18.42 N), located in Lingshui Li Autonomous County, Hainan Province, China. Based on the characteristics of the rice plant, it can be categorized into four main parts: spikes, stalks, leaves, and roots, as illustrated in Figure 1a.

2.2. Biological Characteristics Measurement

2.2.1. Structural Measurement

To simplify the modeling process, the leaves and roots of the rice plant were excluded, resulting in a division of the plant into two main components, the spike and the stalk. The spike, depicted in Figure 1b, is characterized by its panicle shape and comprises rachis, primary branches, two secondary branches, grain stalks, and grains. The primary branches extend from the rachis, and secondary branches may develop on these primary branches. The grains are attached to the branches via grain stalks, with the majority of the grains predominantly located on the exterior of the branches. Furthermore, Figure 1c illustrates that the stalk is segmented into three stalks, labeled primary, secondary, and tertiary, arranged from top to bottom, with each stalk connected by stalk nodes. The overall shape of the stalk tapers from the bottom to the top, exhibiting a gradually increasing diameter between segments. Measurements of the grains, rachis, primary branches, secondary branches, and stalks were conducted using a vernier caliper (model: DHGDW150S, measuring range: 300 mm, graduation value: 0.01 mm, accuracy: ±0.02 mm, Delixi Electric Co., Ltd., Yueqing, China). The measurement process was repeated 15 times to ensure reliability. As shown in Figure 1d, the three-dimensional dimensions of the grains were recorded for length, width, and thickness. As shown in Figure 1e, in addition to measuring the length of the stalks, it was essential to measure both the inner and outer diameters of the hollow stalk structure to accurately determine the wall thickness. Meanwhile, for the rachis, primary branches, and secondary branches, only the length and outer diameter were measured. The structural parameters of each component of the rice plant are presented in Table A1.

2.2.2. Moisture Content and Density Measurement

Due to the significant differences in density and moisture content among various parts of the rice plant, for measurement the rice plant was categorized into three distinct components: grains, branches, and stalks. The branches encompass the rachis, primary branches, and secondary branches, while the stalks consist of the primary, secondary, and tertiary stalks. The moisture content and density of each rice component were evaluated using the drying method and the drainage method. Detailed results and methodologies for these measurements are presented in Appendix Table A1 and Table A2.

2.2.3. Contact Parameter

To enhance the accuracy of the rice discrete element model and optimize its practical application, it was essential to precisely measure the contact parameters between the various components of the rice and the steel and acrylic materials, as well as the rice structure itself. The primary branches, secondary branches, rachis, and stalks of the rice plant were categorized collectively as stalks, and the contact parameters between the grains, stalks, and the steel and acrylic materials were assessed. The grains and stalks were evenly adhered to different cardboards to create grain boards and stalk boards, which were utilized to measure the contact coefficient of the plant structure itself, as illustrated in Figure 2a,b. The coefficients of static friction, rolling friction, and restitution were measured using the grade test, slope rolling test, and drop test, respectively, as illustrated in Figure 2c–e. Detailed measurement methods and results are shown in the attached Appendix Table A2 and Table A3.

2.3. Mechanical Characteristic Test

As illustrated in Figure 3a, a texture analyzer (model: TMS-Pilot, force resolution: 0.01 g, distance accuracy: 0.001 mm, Food Technology Corporation, Sterling, VA, USA) was employed to investigate the mechanical properties of rice. The analysis focused on the grain stalks, primary branches, spikes, and three kinds of stalks undergoing tensile and shear tests. With the exception of the grain stalk, the length of the samples to act as structures in the other parts of the rice plant was 40 mm.

• (1). Tensile test

The tensile fixture comprised upper and lower clamping mechanisms. The opening of the clamp was adjusted by rotating the knob to secure the object being tested. The lower clamping mechanism remained stationary, while the upper clamping mechanism stretched the sample by moving upward. The tensile tests on the rice plant structures were conducted using a texture analyzer equipped with a tensile fixture, as depicted in Figure 3b. The two ends of the sample were secured by a clamping device, while the clamping section, which housed a force sensing element, moved upward at a constant speed of 30 mm/min. Preliminary pre-tests determined the clamping end displacement to be set at 10 mm. As the clamping end ascended, the sample eventually fractured. Each set of tests was repeated five times. The tensile forces are documented in Appendix Table A4, which presents the maximum, minimum, mean, and standard deviation values.

• (2). Shear test

The shear device consisted of a shearing piece and a bottom plate. The shearing piece was oriented perpendicularly to the bottom plate. Additionally, the bottom plate featured guidance grooves. There was a slit in the bottom plate to cut the test sample with the shearing piece. The shearing tests were conducted using a texture analyzer and shearing device, as illustrated in Figure 3c. Prior to initiating the test, it was essential to define the zero level. Preliminary tests indicated that when the shearing device was positioned at −5 mm, the stalk was completely severed. The thickness of the shearing piece was set at 1 mm. The shearing device, along with the force sensor, descends at a loading speed of 60 mm/min until the sample is cut. Each test set was repeated five times. The shear forces are documented in Appendix Table A4, which presents the maximum, minimum, mean, and standard deviation values.

2.4. Establishment of Rice Model

Based on the previously established structural characteristic measurements and mechanical characteristic tests, the whole-rice-plant model was established using the multi-dimensional meta-particle arrangement method. The modeling process is illustrated in Figure 4.

2.4.1. Discrete Element Model of Rice Plant

In the study of rice plant structure, the average values of each structural parameter were selected as the foundation for modeling.

• (1). Establishment of spike simulation model

Prior to modeling, the counts of primary branches on the spike ranged from 7 to 13, while the number of grains on the primary branches varied between 4 and 10. To facilitate the establishment of the model, the number of primary branches on the spike head was established as 10, which included 2 top primary branches and 8 side primary branches. The side primary branches contained 7 grains each, whereas the top primary branches had 6 grains. Therefore, the spike model consists of a rachis, 10 primary branches, and 80 grains and grain stalks, among which the grains and primary branches were connected through the grain stalks, as shown in Figure 5a,b.

The diameter of each structure of the rice spike model particles was determined according to the parameters measured in Appendix Table A1. The average diameter and length of the primary branch were 0.82 mm and 6.43 cm, respectively. The average diameter and length of the rachis were 1.99 mm and 8.46 cm, respectively. The average diameter and length of the grain stalk were 0.59 mm and 0.22 cm, respectively. To facilitate the calculation of particle coordinates, the diameter of the primary branch model particles was set to 1 mm, with a length of 32 mm, comprising a total of 64 single particles. The diameter of the rachis particles was 2 mm, with a length of 84 mm, consisting of 42 particles. For the grain stalk model, it had initially been intended to be composed of 4 particles, each with a diameter of 0.5 mm. However, this configuration would have resulted in contact between grains and other primary branches, leading to unreasonable movement or particle loss. Therefore, to ensure the stability of the rice model, the grain stalk model consisted of a single particle with a diameter of 0.5 mm.

After measuring the dimensions of the grains, the average length, width, and thickness were found to be 8.29 mm, 3.89 mm, and 2.65 mm, respectively. Shuyounuo No. 81 is a variety of japonica rice. Given the short and thick characteristics of japonica rice grains, the grain model was constructed using two particles with a diameter of 2.5 mm, two particles with a diameter of 3.5 mm, and one particle with a diameter of 4 mm. In addition, the primary branches were arranged in a staggered pattern around the rachis, with an angle of 90 degrees between adjacent branches in the radial section. The angle between the primary branch and the main rachis was denoted as α, while the angle between the grain and the branch was represented as β.

• (2). Establishment of stalk simulation model

As shown in Figure 5c,d, rice plants consist of three groups of stalks. To enhance the resemblance of the rice model to the actual rice structure, each stalk features a hollow design. The diameter of the rice stalk model particles was determined according to the parameters measured in Appendix Table A1. The average values of the outer diameter, inner diameter, and length of the primary stalk were 1.96 mm, 1.28 mm, and 24.29 cm, respectively. For the secondary stalk, the average values of the outer diameter, inner diameter, and length were 3.08 mm, 2.03 mm, and 25.12 cm, respectively. The average values of the outer and inner diameters of the tertiary stalk were 4.02 mm, 2.77 mm, and 23.01 cm, respectively. To reduce the complexity of the rice stalk model, the outer diameter, inner diameter, and length of the primary stalk model were set at 2.0 mm, 1.2 mm, and 24.28 cm, respectively. For the secondary stalk model, the corresponding values were set at 3.0 mm, 2.0 mm, and 25.10 cm. The outer diameter, inner diameter, and length of the tertiary stalk model were set at 4.0 mm, 2.8 mm, and 23.04 cm, respectively. Therefore, the primary stalk model consisted of 607 interconnected hollow rigid rings. Each rigid ring was composed of 13 particles, each with a diameter of 0.4 mm, evenly arranged around a circle with a diameter of 1.6 mm. The secondary stalk model comprised 502 interconnected hollow rigid rings, each consisting of 16 particles with a diameter of 0.5 mm, evenly arranged around a circle with a diameter of 2.5 mm. The tertiary stalk model was made up of 384 interconnected hollow rigid rings, where each rigid ring contained 18 spherical particles, each measuring 0.6 mm in diameter, evenly arranged around a circle with a diameter of 3.4 mm.

• (3). Solving grain coordinates of whole-rice-plant model

Using the corresponding model parameters, the particle coordinates for each section of the rice plant model were subsequently calculated. The center of the bottom ring of the tertiary stalk was designated as the origin of the coordinates, thereby establishing a three-axis global coordinate system. Through a series of coordinate translation and rotation transformations, the particle coordinates for each part of the plant model were determined. For a comprehensive outline of the steps involved in calculating the particle coordinates of the rice plant model, please refer to Appendix A.

2.4.2. Model of Particle Bonding

The bonding force between the particles was described by Hertzian contact theory. Notably, there was a significant distinction between the cross-sectional area of the solid particle stalk and the actual cross-sectional area. The cross-sectional area of the hollow ring stalk developed in this study provides a more accurate representation of the normal and tangential characteristics of the real stalk compared to the solid particle stalk. The Hertz–Mindlin band-bonding contact model was employed for interactions between particles, wherein forces and moments were transmitted at the contact points, as illustrated in Figure 6.

This article primarily investigated the particle bonding effects among flexible hollow stalk particles. For a hollow flexible stalk, the parallel bond between the two rings was considered equivalent to a hollow Bernoulli beam. Based on the principles of material mechanics, the normal stiffness and tangential stiffness of the beam per unit area of the ring section were defined.

(1)$k_n={\displaystyle \frac{E_P}{L}}$

(2)$k_{\tau }={\displaystyle \frac{G_P}{L}}$

where k_n is the normal stiffness of the hollow particle, N/m; k_τ is the tangential stiffness of the hollow particle, N/m; E_P is the elastic modulus of the hollow particle, Pa; G_P is the shear modulus of the hollow particle, Pa; and L is the distance between the centers of the two hollow particles, m.

When the particle shape and packing were fixed, the macroscopic Poisson ratio was related to the ratio of k_n to k_τ [36]. Consequently, the normal–shear stiffness ratio of the hollow particles was established to be equal to the actual elastic–shear modulus ratio.

(3)${\displaystyle \frac{k_n}{k_{\tau }}}={\displaystyle \frac{E}{G}}=2\left(1+\mu \right)$

where E is the actual elastic modulus of rice, Pa; G is the true shear modulus of rice, and Pa; and μ are the Poisson ratio of rice.

The overall normal stiffness of the stalks was established through the series connection of parallel bonds between N rings [37]. Consequently, the overall normal stiffness of the stalks can be expressed as follows.

(4)${\displaystyle \frac{1}{k_n^i}}=\sum _{\dot{i}=n}^{N-1}{\displaystyle \frac{1}{A_{\nu }{k}_n}}={\displaystyle \frac{N-1}{A_{\nu }{k}_n}}$

(5)$k_n^i={\displaystyle \frac{A_{\nu }{k}_n}{N-1}}$

The overall tangential stiffness of the stalk can be obtained by connecting Equations (3) and (5).

(6)$k_{\tau }^i={\displaystyle \frac{A_{\nu }{k}_n}{2\left(1+\mu \right)\left(N-1\right)}}$

where k_nⁱ is the normal stiffness of the whole-stalk model, N/m; k_τⁱ is the tangential stiffness of the whole-stalk model, N/m; A_v is the cross-sectional area of the hollow stalk model, m²; and N is the number of hollow stalk particles.

After conducting a tensile test using a texture analyzer, the slope of the tensile curve represented the normal stiffness at that specific length. The normal stiffness for each structure of the rice plant is detailed in Appendix Table A4. It was assumed that the normal stiffness of the rice plant structure in the tensile test sample was equivalent to the normal stiffness of the stem particle model of the same length. We can substitute the actual normal stiffness of the structure into Equation (5). This allows us to determine the normal stiffness per unit area of the parallel bonds for each section of the rice structure, as illustrated in Equation (7).

(7)$k_n^r=k_n^i={\displaystyle \frac{A_v{k}_n}{N-1}}$

Transform Equation (7) to obtain Equation (8).

(8)$k_n={\displaystyle \frac{\left(N-1\right)k_n^r}{A_{\nu }}}$

where k_n^r is the actual normal stiffness of the stalk, N/m.

When the bond parameters between particles reach critical values, the bonding connections break, causing the particles to revert to their discrete characteristics. The maximum tensile and shear forces for each structure of the rice plant can be determined through tensile and shear tests, with detailed results recorded in Appendix Table A5 and Table A6, respectively. The critical normal stress and critical shear stress can be calculated using the established Equations (9) and (10).

(9)${\sigma }_{max}={\displaystyle \frac{F_n^{max}}{A_r}}$

(10)${\tau }_{max}={\displaystyle \frac{F_{\tau }^{max}}{A_r}}$

where σ_max is the critical tensile stress, N/m²; F_n^max is the maximum tensile force, N; F_τ^max is the maximum shear stress; and A_r is the actual cross-sectional area of the stalk, m².

Using the aforementioned formula, the normal stiffness per unit area, the shear stiffness per unit area, the critical tensile stress, and the critical shear stress for each structure of the rice plant can be determined. The simulated bonding parameter ranges of each structure of the rice plant model are shown in Appendix Table A7.

2.5. Experimental Design for Optimization of Bending Characteristics of Hollow Stalk Model

2.5.1. Three-Point Bending Test Scheme

To calibrate the bending characteristics, the three-point bending test was conducted on the rice stalks. It was observed that the bending resistance strength of the rice stalk nodes was greater than that of the stalk, and the bending characteristics of the stalk nodes were less pronounced. Consequently, the bonding parameters at the stalk nodes were set to a mean value, and this study focuses exclusively on the bending characteristics of the stalk. Samples measuring 150 mm in length were selected from the middle of the primary, secondary, and tertiary stalks, respectively, with a span of 130 mm between the support points at both ends, as illustrated in Figure 7a. Prior to the test, pre-tests indicated that the stalk would undergo plastic bending if the descent depth exceeded 20 mm. Therefore, to optimize the bending characteristics of the stalk under radial loading while avoiding compression damage in the radial direction, the descent depth was limited to 15 mm. During the test, the initial position of the bending probe was set 10 mm above the support surface of the sample, and it was moved vertically downward toward the center of the sample. The loading speed was maintained at 60 mm/min, causing the stalks to bend under the applied load. Force changes began to be recorded upon contact between the probe and the sample, and when the probe descended to a depth of 15 mm, it was retracted at a speed of 180 mm/min, concluding the test. Each experimental group was repeated five times, and the average value was calculated. The curve of actual bending resistance force is shown in Figure 7c.

As illustrated in Figure 7b, the three-point bending test was employed to simulate the deflection of rice stalks using EDEM (2020 version, Altair Engineering Inc., Troy, Michigan, USA). The bending system, constructed with the same device model, was prepared in SolidWorks (version 2021, Dassault Systemes Inc., Vélizy-Villacoublay, France). During the simulation process, the material properties, the spacing between support ends, the probe loading speed, and the positional relationships were configured to align with the actual bending test conditions. With the exception of five bonding parameters, the remaining simulation parameters were averaged, while the unmeasured simulation parameters were established by consulting the literature [38,39,40,41,42], as detailed in Table 1. The simulation duration was set to 35 s, with a time step of 2 × 10⁻⁶, and the test data were recorded every 0.01 s. The curve of the simulated bending resistance force is shown in Figure 7d.

2.5.2. Bending Characteristics Optimization Test Scheme

This study employed the Plackett–Burman screening test, the steepest ascent test, and the Box–Behnken optimization test to optimize the bending characteristics for the three kinds of stalk. The five bonding parameters—normal stiffness per unit area, shear stiffness per unit area, critical normal stress, critical shear stress, and bonded disk radius—were utilized as test factors, while the bending resistance force of the stalks served as the test index for parameter optimization. The range of each experimental factor was established based on prior experiments and a literature review [43,44], and was coded using the +1 and −1 levels to represent the maximum and minimum values, respectively, as detailed in Table 2.

2.6. Experimental Design for Optimizing Threshing Characteristics of Rice Model

2.6.1. Impact Threshing Test

To more accurately reflect the deformation and threshing conditions of rice plants under impact loads, we optimized the threshing characteristics of the rice model through impact threshing tests. As illustrated in Figure 8a, the experimental apparatus consists of a pendulum, base plate, support device, and fixed installation. The pendulum weighs 300 g and has a length of 250 mm. To determine the appropriate pendulum angle, a preliminary experiment was conducted prior to the formal testing. The rice plant was placed in the support device, and the initial angle of the pendulum was adjusted accordingly. During the experiment, the pendulum angle was set to 15°, 30°, 45°, 60°, 75°, and 90° to perform the impact tests. Upon release, the pendulum collides with the rice spike, resulting in the elastic bending of the stalk and the subsequent detachment of the grains. The number of grains that fell off after the impact was recorded, and the threshing rate was calculated. Each test at a given angle was repeated ten times, and the average value was computed. The formula for calculating the threshing rate is as follows:

(11)$P={\displaystyle \frac{N_1}{N_2}}\times 100\%$

where N₁ is the number of fallen grains, and N₂ is the total number of grains.

The threshing rates at various pendulum angles are presented in Appendix Table A8. Preliminary experiments indicated that as the pendulum angle increased, the number of grains dislodged by the impact also increased. This phenomenon occurred because a higher lifting angle resulted in greater gravitational potential energy for the pendulum. During the impact process, this gravitational potential energy was converted into increased kinetic energy, thereby enhancing the impact effect. To facilitate statistical analysis and observation, a pendulum angle of 45° was selected. The impact threshing test ignores the effects of air resistance and friction.

As illustrated in Figure 8b, the EDEM–Recurdyn coupling method was employed to perform the impact simulation test. Initially, in Recurdyn (2023 version, FunctionBay Inc., Gyeonggi-do, KOR), the angle of the pendulum was set to 45° and allowed to swing downward under the influence of gravity. Subsequently, the rice plant model was created in EDEM, configured to generate two rice plant models within support device over a duration of 0.1 s. The pendulum, influenced by gravity, descended and struck the rice plant, resulting in the grains detaching. The simulation time was 0.3 s, with a time step of 2 × 10⁻⁶, and test data were recorded every 0.01 s. Each group of tests was conducted ten times. The number of grains that fell off after the impact was recorded, the threshing rate of the grains was calculated, and the average value was determined.

2.6.2. Threshing Characteristics Optimization Test Scheme

The Plackett–Burman screening test, steepest ascent test, and Box–Behnken optimization test were employed to optimize the bonding mechanical parameters of the grain stalk. Five bonding parameters—normal stiffness per unit area, shear stiffness per unit area, critical tensile stress, critical shear stress, and bonded disk radius—were selected as experimental factors. These parameters were optimized based on the experimental index of the threshing rate of the rice model under impact load. To minimize the number of tests and experimental factors, the bonding parameters for grain stalk–grain, grain stalk–primary branch, and grain stalk–grain stalk were considered equal, while the bonding parameters for grain stalk–grain stalk were optimized. The simulation parameters for other structures were averaged. The range of test factors was established through preliminary experiments, as detailed in Table 3.

2.7. Design of the Rice Model Stack Angle Verification Test

The stack angle test serves to evaluate the stacking condition of the rice plant under gravitational forces and to verify the contact characteristics of the rice plant model. The verification tests were conducted to assess the stacking angles between grain–grain, stem–stem, and grain–stem. The secondary stalk was selected for the test.

2.7.1. Grain–Grain Stack Angle Test

A hollow acrylic tube with a diameter of 30 mm, a height of 100 mm, and a wall thickness of 3 mm was placed vertically on a flat surface. Two thousand grains were introduced into the tube, and the tube was gradually lifted. The grains rolled down due to the force of gravity until they were no longer in contact with the tube, at which point their movement ceased. The test results are presented in Figure 9a. This test was conducted for 3 s, and the average value was calculated. In the simulation test, as depicted in Figure 9d, the particle factory was configured at a height of 50 mm, with a particle generation rate of 4000 particles per second. The particles accumulated in the acrylic tube following their generation. After 1 s, the acrylic tube was slowly lifted at a speed of 0.1 m/s. The total simulation duration was 5 s, with a time step of 2 × 10⁻⁶, and test data were recorded every 0.01 s. For the stack angle analysis, Matlab (2020b version, MathWorks Inc., Natick, Massachusetts, USA) was utilized to process the stacked images sequentially through grayscale conversion, binarization, and edge extraction. The extracted image information from the RC coordinate system was transformed into the XY coordinate system to derive the stacked unilateral contour line, which was subsequently performed with the linear fitting curve. The grain–grain accumulation images were processed through grayscale conversion and binarization, as shown in Figure 9b,e, with the linear fitting curves illustrated in Figure 9c,f.

2.7.2. Stalk–Stalk Stack Angle Test

In the stalk–stalk stack angle test, the stalks were processed into 15 mm segments, resulting in a total of 200 segments. These stalks were placed in an acrylic box with dimensions of 80 mm × 30 mm × 60 mm. One side wall of the box was gradually raised to allow the stalks to move under the influence of gravity. The test results are presented in Figure 10a. For the simulation test, a particle factory was established at a height of 60 mm from the bottom of the container. Particles accumulated in the container due to gravitational forces, with a stalk particle generation rate of 1000 particles per second. After 1 s, the acrylic plate was slowly lifted at a speed of 0.1 m/s and subsequently allowed to fall, causing the particles to accumulate under gravity. The simulation test results are illustrated in Figure 10d. The total simulation time was 5 s, with a time step of 2 × 10⁻⁶, and test data were recorded every 0.01 s. The stalk–stalk stack image was processed using grayscale conversion and binarization, as shown in Figure 10b,e, while the linear fitting curves are depicted in Figure 10c,f.

2.7.3. Grain–Stalk Stack Angle Test

For the grain–stalk stack angle test, 200 sections of 15 mm length stalks and 1000 grains were placed into an acrylic box with dimensions of 80 mm × 30 mm × 60 mm. The testing method was consistent with the stalk–stalk test. The results are presented in Figure 11a,d. The grain–stalk stack images were processed using grayscale conversion and binarization, as illustrated in Figure 11b,e, while the linear fitting curves are depicted in Figure 11c,f.

3. Results and Discussion

3.1. The Stalk Bending Characteristics Calibration

3.1.1. Analysis of the PB Test Plan and the Steepest Ascent Test for the Stalk Bending Resistance Force

Table 4 presents the Plackett–Burman test plan and results for calibrating the stem bending characteristic parameters of three rice groups. Regression analysis on the test results was performed using Design-Expert (version 13, Stat-Ease, Inc., Minneapolis, Minnesota, USA), with the findings detailed in Table 5. The p-values of the regression model were all below 0.05, indicating that the regression equation, which utilizes bending resistance force as the test index, accurately reflects the actual conditions.

The regression model of the primary stalk bending resistance force y₁ was as follows:

(12)$y_1=1.02+0.0022x_1+0.0912x_2+0.0362x_3+0.0822x_4+0.0685x_5$

The regression equation indicates that the factors influencing the bending resistance force at the primary stalk, in order of significance, were x₂, x₄, x₅, x₃, and x₁. To assess their impact on bending resistance force, steepest ascent tests were conducted on x₂, x₄, and x₅, while x₁ and x₃ were held constant at their average values during the testing process.

The regression model of the secondary stalk bending resistance force y₂ was as follows:

(13)$y_2=2.33367+0.127x_1-0.045x_2+0.05333x_3+0.10133x_4+0.087667x_5$

The regression equation indicated that the factors influencing the bending resistance force at the secondary stalk, in order of significance, were x₅, x₁, x₄, x₃, and x₂. To assess their impact on bending resistance force, steepest ascent tests were conducted on x₁, x₄, and x₅, while x₂ and x₃ were averaged during the testing process.

The regression model of the tertiary stalk bending resistance force y₃ was as follows:

(14)$y_3=4.759-0.08633x_1+0.20933x_2+0.181337x_3+0.024x_4+0.102x_5$

According to the regression equation and significance analysis, the factors influencing the bending resistance force at the tertiary stalk, in order of impact, were x₂, x₃, x₅, x₁, and x₄. To assess their effects on bending resistance force, steepest ascent tests were conducted on x₂, x₃, and x₅, while _x1 and x₄ were averaged during the testing process.

The results of the simulation steepest ascent test conducted in the stalk are presented in Table 6. For the primary stalk steepest ascent test, the simulation results for the second group (1.097 N) were in close proximity to the target value obtained from the actual test (1.013 N). The parameter values of the second set of simulations should be designated as the center of mass, while the first and third sets were assigned a low level (−1) and high level (+1), respectively.

For the secondary stalk steepest ascent test, the simulation results for the third group (2.576 N) were quite close to the target value obtained from the actual test (2.464 N). Consequently, designating the parameter values from the third set as the centroid, while the values from the second and fourth sets were assigned a low level (−1) and high level (+1).

For the tertiary stalk steepest ascent test, the simulation results for the third group (4.705 N) were notably close to the target value obtained from the actual test (4.593 N). The third set of simulation parameter values were used as the centroid, while the second and fourth sets of values were designated a low level (−1) and high level (+1), respectively.

The level coding table for the test factors of the three stalk bending tests is shown in Table 7, and the results of the Box–Behnken test are shown in Table 8.

3.1.2. Analysis of the Box–Behnken Test for the Stalk Bending Resistance Force

The variance analysis of the experimental results, as shown in Table 9, indicated that the p-values of the regression models were all less than 0.01, while the p-values of the residual models were all greater than 0.05. This suggested that the regression equations, with bending resistance force as the experimental index, exhibited good fitting performance. As presented in Figure 12, the correlation coefficients (R) for the experiments were all greater than 0.98, indicating a strong linear relationship between the experimental factors and the experimental index. Additionally, the coefficients of determination (R²) were all greater than 0.97, demonstrating a high degree of fit for the regression models. Furthermore, the adjusted coefficients of determination (Adj R²) were all greater than 0.935, which implied that the experimental factors in the models accounted for more than 93.50% of the variation in the experimental index.

In the primary stalk bending test, x₂, x₄, x₅, and x₅² had a highly significant effect on the bending resistance force (p < 0.01). x₄x₅ and x₂² had a significant effect on the bending resistance force (0.01 < p < 0.05). x₂x₄, x₂x₅, and x₄² had no significant effect on the bending resistance force (p > 0.05). The order of influence of the experimental factors on the primary stalk bending resistance force, from largest to smallest, was as follows: shear stiffness per unit area, critical shear stress, and bonded disk radius. After removing the insignificant terms, the actual regression equation for the primary stalk bending resistance force was obtained, as shown in Equation (15).

(15)$\begin{array}{c}{y}_1=-5.12252+8.91703\times {10}^{-12}{z}_2-0.059895z_4-6.12928z_5\\ +0.387368z_4{z}_5-2.99284\times {10}^{-24}{z_2}^2-9.20084{z_5}^2\end{array}$

The response surface plots illustrating the effects of shear stiffness per unit area, critical shear stress, and bonded disk radius on the bending resistance force of the primary stalk are shown in Figure 13. It can be observed that, as the shear stiffness per unit area and critical shear stress increased, the bending resistance force of the stalk model gradually increased. However, as the bonded disk radius increased, the bending resistance force of the stalk model first increased and then decreased. By utilizing the optimization module in Design-Expert and setting a target bending force of 1.013 N (obtained from the calibration experiments) as the objective, the parameters of the three factors were optimized. The optimal fitting parameters obtained were a shear stiffness per unit area of 1.3 × 10¹² N/m³, a critical shear stress of 36.95 MPa, and a bonded disk radius of 0.437 mm. To validate the optimization results, three simulation tests were conducted using the optimized significant parameter values and the average values of the insignificant parameters. The bending resistance forces obtained from the simulations were 1.086 N, 1.061 N, and 1.034 N, with an average value of 1.0603 N. The relative error compared to the actual measured value was 4.46%.

In the secondary stalk bending test, x₁, x₄, x₅, x₁², and x₅² had a highly significant effect on the bending resistance force (p < 0.01). x₄x₅ had a significant effect on the bending resistance force (0.01 < p < 0.05). x₁x₄, x₁x₅, and x₄² had no significant effect on the bending resistance force (p > 0.05). The order of influence of the experimental factors on the secondary stalk bending resistance force, from largest to smallest, was as follows: normal stiffness per unit area, critical shear stress, and bonded disk radius. After removing the insignificant terms, the actual regression equation for the secondary stalk bending resistance force was obtained, as shown in Equation (16).

(16)$\begin{array}{c}{y}_2=-33.47524+6.95909\times {10}^{-11}{z}_1-0.151767z_4-2.80166z_5\\ +0.348279z_4{z}_5-3.21444\times {10}^{-23}{z_1}^2-3.7985{z_5}^2\end{array}$

The response surface plots illustrating the effects of normal stiffness per unit area, critical shear stress, and bonded disk radius on the bending resistance force of the secondary stalk are shown in Figure 14. It was observed that as the normal stiffness per unit area and critical shear stress increased, the bending resistance force of the secondary stalk model gradually increased. Meanwhile, as the bonded disk radius increased, the bending resistance force first increased and then decreased. The optimization module of Design-Expert was employed to optimize the parameters of the three factors, with the target set to a bending resistance force of 2.464 N, obtained from calibration tests. The optimal parameter values were determined as follows: normal stiffness per unit area of 1.09 × 10¹² N/m³, critical shear stress of 23.36 MPa, and bonded disk radius of 0.6 mm. To validate the optimization results, three simulation experiments were conducted using the optimized significant parameter values and the average values of the non-significant parameters. The bending resistance forces obtained from the simulations were 2.513 N, 2.429 N, and 2.158 N, with an average value of 2.3667 N. The relative error compared to the actual measured value was 3.95%.

In the tertiary stalk bending test, x₂, x₃, x₅, and x₅² had a highly significant effect on the bending resistance force (p < 0.01). x₂x₃ and x₂² had a significant effect on the bending resistance force (0.01 < p < 0.05). x₂x₅, x₃x₅, and x₃² had no significant effect on the bending resistance force (p > 0.05). The order of influence of the experimental factors on the tertiary stalk bending resistance force, from largest to smallest, was as follows: shear stiffness per unit area, critical normal stress, and bonded disk radius. After removing the insignificant terms, the actual regression equation for the tertiary stalk bending resistance force was obtained, as shown in Equation (17).

(17)$\begin{array}{c}{y}_3=16.16817-1.11429\times {10}^{-12}{z}_2-1.57661z_3+4.542998z_5\\ +6.80952\times {10}^{-12}{z}_2{z}_3-2.51429\times {10}^{-22}{z_2}^2-2.87663{z_5}^2\end{array}$

The response surface plots illustrating the effects of shear stiffness per unit area, critical normal stress, and bonded disk radius on the bending resistance force of the secondary stalk are shown in Figure 15. As the shear stiffness per unit area and critical normal stress increased, the bending resistance force of the stalk model gradually increased. However, with the increase in the bonded disk radius, the bending resistance force first increased and then decreased. Using the optimization module in the Design-Expert software, the parameters of the three factors were optimized with the target set to a bending resistance force of 4.593 N, obtained from the calibration tests. The optimal parameter values were determined as follows: a shear stiffness per unit area of 2.6 × 10¹¹ N/m³, a critical normal stress of 21.73 MPa, and a bonded disk radius of 1.09 mm. To validate the optimization results, three simulation experiments were conducted using the optimized significant parameter values and the average values of the non-significant parameters. The bending resistance forces obtained from the simulations were 4.260 N, 4.661 N, and 4.512 N, with an average value of 4.4777 N. The relative error compared to the actual measured value was 2.51%.

3.2. Analysis of Threshing Characteristic Calibration Test

3.2.1. Analysis of PB Test Scheme and Steepest Climb Test for Model Threshing Rate

The Plackett–Burman experimental design and results for the calibration of threshing characteristic parameters of the entire rice model are presented in Table 10. The regression analysis of the experimental results is shown in Table 11. The p-value of the regression model was less than 0.05, indicating that the regression equation fitted to the experimental data for the threshing rate aligned well with the actual observations.

The regression model of the threshing rate y₃ was as follows:

(18)$y_4=38.78-1.13x_1-1.647x_2-1.83x_3-0.4967x_4-1.33x_5$

According to the regression equation and significance analysis, the factors influencing the threshing rate of the model, in descending order of importance, were x₃, x₂, x₅, x₁, and x₄. To further determine their effects on the threshing rate, a steepest ascent experiment was conducted for x₁, x₂, x₃, and x₅. During the experiment, x₄ was held at its average value.

The results of the steepest ascent experiment, using the model threshing rate as the experimental indicator, are presented in Table 12. Among the results, the simulation outcome of the third group (38.14%) was the closest to the target value from the actual experiment (38.71%). Therefore, the parameter values of the third group were designated as the centroid, and parameter optimization experiments were conducted within the numerical range defined by the second and fourth groups. The level coding table for the experimental factors is shown in Table 13, and the design and results of the Box–Behnken experiment are presented in Table 14.

3.2.2. Analysis of the Box–Behnken Test for Model Threshing Rate

The ANOVA for the model threshing rate results is shown in Table 15. The p-value of the regression model was less than 0.01, while the p-value of the residual model was greater than 0.05, which indicated that the regression equation fitted well to the experimental data, with the threshing rate of the rice model as the experimental indicator. The correlation coefficient (R) of the experiment was 0.977, suggesting that there was a strong linear relationship between the experimental factors and the experimental indicator. The coefficient of determination (R²) was 0.9548, which demonstrated that the regression model had a good fit. Furthermore, the adjusted coefficient of determination (Adj R²) was 0.9153, which indicated that the experimental factors in the model explained more than 91.53% of the variation in the experimental indicator.

In the model threshing rate experiment, x₁, x₂, x₃, x₅, x₁x₂, x₁x₃, x₂x₃, x₁², and x₅² had extremely significant effects on the threshing rate of the model (p < 0.01). x₂x₅ had a significant effect on the threshing rate (0.01 < p < 0.05). Additionally, x₁x₅, x₃x₅, x₂², and x₃² showed no significant effects on the model threshing rate (p < 0.05).

The order of influence of the experimental factors on the threshing rate of the model, from greatest to least, was as follows: critical tensile stress, shear stiffness per unit area, bonded disk radius, and normal stiffness per unit area. After removing the insignificant terms, the actual regression equation for the threshing rate of the model, under a swing angle of 45°, was obtained as shown in Equation (19).

(19)$\begin{array}{c}{y}_4=-215.80885-2.6652\times {10}^{-8}{z}_1-2.39913\times {10}^{-8}{z}_2-18.17089z_3+11.32429z_5\\ +1.89412\times {10}^{-18}{z}_1{z}_2+1.05385\times {10}^{-9}{z}_1{z}_3+2.67873\times {10}^{-9}{z}_2{z}_3\\ -1.95765\times {10}^{-8}{z}_2{z}_5+9.63942\times {10}^{-19}{z_1}^2-1.00525{z_5}^2\end{array}$

Assuming that x₃ = x₅ = 0, x₂ = x₅ = 0, x₁ = x₅ = 0, and x₁ = x₃ = 0, we can derive the influence rules of the interaction factors x₁x₂, x₁x₃, x₂x₃, and x₂x₅ on the model threshing rate y₄. The results are illustrated in Figure 16.

With the increase in normal stiffness per unit area, the threshing rate of the model first increased and then decreased. This phenomenon occurred because the increase in normal stiffness per unit area enhanced the bonding force between particles, thereby increasing the resistance to impact of the particles. However, as the normal stiffness per unit area continued to increase, the bonding stiffness between particles also increased, which made the bonding brittle. This brittleness caused the bonding to break more easily under stress.

As the shear stiffness per unit area increased, the threshing rate of the model gradually decreased. This was because the increase in shear stiffness improved the torsional resistance between particles, partially counteracting the impact forces during collisions.

When the critical tensile stress increased, the threshing rate decreased. This was because the increase in critical tensile stress raised the ultimate breaking limit of the bonding between particles, allowing the particles to withstand greater impact forces without breaking.

As the bonded disk radius increased, the threshing rate decreased. This was because the increase in the bonded radius expanded the range of the bonding force, which enhanced the bonding force to some extent. This increase in bonding force also improved the toughness of the bonding between particles, making the bonds harder to break.

With the target set to the calibrated experimental average impact threshing rate of 38.71%, parameter optimization was performed. The obtained optimal fitting parameters were a normal stiffness per unit area of 7.04 × 10⁹ N/m³, a shear stiffness per unit area of 2.91 × 10⁹ N/m³, a critical tensile stress of 4.77 MPa, and a bonded disk radius of 2.88 mm. To validate the optimized parameters, three simulation verification experiments were conducted using the optimized significant parameter values and the average values of the non-significant parameters. The resulting impact threshing rates of the model were 37.11%, 38.84%, and 38.02%, with an average value of 37.99%. The relative error compared to the actual measured value was 1.86%, indicating that the model could effectively simulate the grain detachment behavior.

3.3. Analysis of Stack Angle Verification Test Results

From the side view of the stack angle verification test, the sample being tested slides downward due to the force of gravity and accumulates. In the front view, the sample slides forward, losing the constraints imposed by the baffles on both sides, and spreads outward in a fan-shaped manner. The results of the simulation tests align with the actual experimental observations. The comparison of stack angle verification tests involving grain–grain, stalk–stalk, and grain––stalk interactions is illustrated in Figure 17, revealing relative errors of 4.11%, 7.52%, and 5.73%, respectively. Notably, the stalk–stalk relative error was the largest, attributed to the significant impact of moisture content on the test results. The experiments indicate a positive correlation between moisture content and stack angle; as the moisture content of rice increases, the friction characteristics of the stalk surface are enhanced, resulting in a larger stacking angle. Furthermore, the stack angles for stalk–stalk and stalk–grain interactions were slightly larger than those observed in the simulation tests, which may be related to the structural characteristics of the stalk model. The stalk model consists of many spherical particles, rather than being a continuous and smooth entity, exhibiting pronounced regular undulations and distinct concave and convex features. In the grain–grain contact test, the actual test results were slightly larger than those of the simulation test due to the rough surface of the actual grains, which feature small burr structures at their ends, thereby increasing their friction properties.

3.4. Discussion

In current rice plant modeling methods, rice plant models typically consist of two components: a grain model and a stalk model [20,21,22,23]. However, this model lacks overall coherence and is inadequate for simulating continuous deformation processes. Furthermore, existing entire rice plant models significantly deviate from the actual rice plant [32,33,34,35]. The particle arrangement method employed in rice plant models fails to adequately support the simulation and analysis of processes such as the separation of grains and stalks, resulting in a lack of simulation credibility. Consequently, there remained a need to further enhance the completeness and accuracy of rice models. By analyzing the characteristics of the rice plant and measuring key physical parameters, a segmented combined whole-rice-plant modeling method was proposed, establishing a whole-rice-plant model with flexibility and threshing characteristics. The feasibility and accuracy of the model were validated through five groups of tests and simulation results. The analysis of the results from the three stalk bending simulation tests revealed that the accuracy errors for the primary, secondary, and tertiary stalks were 4.46%, 3.95%, and 2.51%, respectively. Shi et al. established a hollow stalk model, and the error rate of the model in a bending test was 10.21% [17]. Su et al. developed a flexible solid stalk model and verified it by a bending test where the error rate of the model was 6.10% [19]. The maximum error rate in this study was decreased by 56.32% and 26.89% compared to the previous two studies, indicating that the hollow stalk structure model and the ring parallel bonding mechanics parameters employed in this research more accurately reflect the curved characteristics of a real stalk. The results of the impact threshing test indicated that, as the pendulum angle increases, the threshing rate of the rice also gradually rises, which aligned with the findings of Sun [44]. Furthermore, at a pendulum angle of 45 degrees, the threshing rate of the model was 1.86%, which effectively demonstrated the feasibility of the rice model concerning its threshing characteristics. The results of the verification test for the stack angle indicated that the error rates for the grain–grain, stalk–stalk, and stalk–grain accumulation angles were 4.11%, 7.52%, and 5.73%, respectively. To enhance the accuracy of rice plant models, we examined the errors produced during simulation tests. The primary source of testing errors was the moisture content, followed by random errors arising during the testing process, which aligns with the findings of Ashtiani [45]. During the tests, the samples were divided into smaller pieces, resulting in significant water loss to the air and considerable fluctuations in moisture content, which affected the test results. Although deviations occur in the simulation tests, their impact on the results was minimal and can be considered negligible in practical applications. The whole-rice-plant model effectively reflected the variations in the characteristics of different parts of the rice plant. It can be utilized to simulate interactions between rice plants and mechanical devices, particularly in terms of the breaking and separation of the rice structure. This provided theoretical support for the study of rice threshing mechanisms during the threshing process.

However, there were still some limitations in the modeling method proposed in this paper. Firstly, while the accurate rice plant model enhances the credibility of harvest simulations, it also increases the simulation duration. Therefore, it was essential to find a balance between simulation accuracy and time efficiency. For the stalk model, the use of single rigid ring particles can be modified to multiple particles to enhance simulation credibility while simultaneously reducing simulation time. Secondly, the modeling process overlooked the impact of the blade and root system, resulting in an incomplete rice model. In full-feed threshing, the rice is cut by the header while the root system remains in the soil. The portion of the rice plant above the root system is fed into the threshing chamber for threshing. The impact of the blade and root system on this process can be disregarded [46]. However, for semi-feed threshing, the model cannot accurately simulate soil–rice–machinery interactions without leaf and root models. In the future, reverse engineering technology can be used to model a rice blade and root system to improve the whole-rice-plant model. Finally, it was essential to develop a rice model that can adjust simulation parameters based on varying moisture content to ensure the general applicability of the rice model.

4. Conclusions

Aiming at solving the problems that rice plant modeling cannot accurately reflect, the difference in rice parts, the characteristics of threshing, and the lack of overall coherence, this paper proposed a method of whole-rice-plant modeling based on segmented hollow stalks. The feasibility and effectiveness of the proposed method were verified by the analysis of experimental and simulation results. The specific research results and discovery were as follows:

1. Based on the characteristics of rice plants, the biological and mechanical characteristics of each segment of a rice plant were studied. On this basis, a segmented hollow-stalk whole-rice-plant model was established by using the multi-dimensional particle arrangement method, and the Hertz–Mindlin contact model was used for the contact mechanics of grain interaction.

2. The stalk bending characteristics and rice threshing characteristics of the established whole-rice-plant model were calibrated. The calibration errors for the three kinds of stalk model in the three-point bending simulation tests were 4.46%, 3.95%, and 2.51%, respectively. In the impact threshing simulation test, the calibration error for the threshing rate was 1.86%, indicating that the model can accurately simulate the bending and threshing behavior of rice plants.

3. The stack angle verification test was conducted on the contact parameters of the rice model. The accuracy errors for the grain–grain, stalk–stalk, and grain–stalk stacking angles were found to be 4.11%, 7.52%, and 5.73%, respectively. Moisture content was the main factor that affects model accuracy and produces experimental errors.

4. The modeling method proposed in this study offers a feasible and effective flexible whole-rice-plant model for simulating rice threshing and cleaning processes. It can be used to simulate the interaction between rice plants and mechanical devices, so as to simulate the threshing and destruction process of rice, and provide a new understanding for the study of rice plant microdynamics. Furthermore, this modeling and parameter calibration method can provide a reference for the development of discrete element models of plant structure for other crop species.

Author Contributions

Conceptualization, R.Y.; data curation, P.W.; formal analysis, P.W., D.C., W.S. and K.X.; funding acquisition, R.Y. and Y.Q.; investigation, Y.Q.; methodology, R.Y., D.C. and W.S.; project administration, R.Y. and Y.Q.; resources, L.C.; software, P.W., L.C. and K.X.; supervision, D.C.; validation, P.W. and W.S.; visualization, Y.Q. and L.C.; writing—original draft, P.W.; writing—review and editing, Y.Q. All authors have read and agreed to the published version of the manuscript.

Institutional Review Board Statement

Not applicable.

Data Availability Statement

The data presented in this study are available on request from the second author at ([email protected]).

Acknowledgments

The authors are pleased to acknowledge the rice farmers for providing rice samples.

Conflicts of Interest

The authors declare no competing interests.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Structural characteristics of rice plants. (a) Whole structure of rice plant. (b) Spike structure. (c) Stalk structure. (d) Grain size. (e) Stalk cross-section.

Figure 2. (a) Grain plate. (b) Stalk plate. (c) Grade test. (d) Slope rolling test. (e) Drop test.

Figure 3. Tests on mechanical properties of rice. (a) Mechanical properties test based on texture analyzer. (b) Tensile test. (c) Shear test.

Figure 4. Modeling process of whole-rice-plant discrete element model.

Figure 5. Discrete element model of different parts of rice plant based on segmentation. (a) Front view of panicle head discrete element model. (b) Side view of panicle head discrete element model. (c) Stalk discrete element model. (d) Discrete element model of whole rice plant.

Figure 6. An illustration of the particle bonding. Where r1 is the physical radius of the particle, mm; r2 is the bonding radius, mm; L is the distance between the two particles, mm; Ft and Fn are the normal force and tangential force of the bonding bond, N; and Mt and Mn are the normal and tangential moments of the bonding bond, N·m.

Figure 7. Three-oint bending simulation test. (a) Actual test. (b) Simulation test. (c) Bending resistance force curve of actual stalk bending test. (d) Bending resistance force curve of simulation stalk bending test.

Figure 8. Impact threshing test. (a) Actual test. (b) Simulation test.

Figure 9. Grain–grain stack angle test. (a) Actual test. (b) Binarization of actual stack images. (c) Actual stack contour linear fitting. (d) Simulation test. (e) Binarization of simulation stack images. (f) Simulation stack contour linear fitting.

Figure 10. Stalk–stalk stack angle test. (a) Actual test. (b) Binarization of actual stack images. (c) Actual stack contour linear fitting. (d) Simulation test. (e) Binarization of simulation stack images. (f) Simulation stack contour linear fitting.

Figure 11. Grain–stalk stack angle test. (a) Actual test. (b) Binarization of actual stack images. (c) Actual stack contour linear fitting. (d) Simulation test. (e) Binarization of simulation stack images. (f) Simulation stack contour linear fitting.

Figure 12. Regression coefficient of the experiment result.

Figure 13. The influence of test factors on the bending resistance force of the primary stalk. (a) The influence of shear stiffness per unit area and critical shear stress on the bending resistance force. (b) The influence of shear stiffness per unit area and bonded disk radius on the bending resistance force. (c) The influence of critical shear stress and bonded disk radius on the bending resistance force.

Figure 14. The influence of test factors on the bending resistance force of the secondary stalk. (a) The influence of normal stiffness per unit area and critical shear stress on the bending resistance force. (b) The influence of normal stiffness per unit area and bonded disk radius on the bending resistance force. (c) The influence of critical shear stress and bonded disk radius on the bending resistance force.

Figure 15. The influence of test factors on the bending resistance force of the tertiary stalk. (a) The influence of shear stiffness per unit area and critical normal stress on the bending resistance force. (b) The influence of shear stiffness per unit area and bonded disk radius on the bending resistance force. (c) The influence of critical normal stress and bonded disk radius on the bending resistance force.

Figure 16. The influence of test factors on the bending resistance force of the model threshing rate. (a) The influence of normal stiffness per unit area and shear stiffness per unit area on the model threshing rate. (b) The influence of normal stiffness per unit area and critical normal stress on the model threshing rate. (c) The influence of shear stiffness per unit area and critical normal stress on the model threshing rate. (d) The influence of shear stiffness per unit area and bonded disk radius on the model threshing rate.

Figure 17. Stack Angle Verification Test Results.

Table 1

Simulation parameter.

Parameters		Value	Source
Modulus of elasticity	Grain	287.35 MPa	Reference [38]
	Stalk	2.8 MPa	Reference [39]
	Steel	2.06 × 105 MPa	Reference [40]
	Acrylic	2.7 GPa	Reference [41]
Poisson ratio	Grain	0.34	Reference [38]
	Stalk	0.4	Reference [39]
	Steel	0.3	Reference [40]
	Acrylic	0.35	Reference [42]
Density	Steel	7800 kg/m3	Reference [40]
	Acrylic	1070 kg/m3	Reference [41]

Table 2

Factors and levels of the PB test for the stalk bonding parameters.

Code	Simulation Parameters	Primary Stalk		Secondary Stalk		Tertiary Stalk
		Low Level (−1)	High Level (+1)	Low Level (−1)	High Level (+1)	Low Level (−1)	High Level (+1)
x1	Normal stiffness per unit area z1 (N/m3)	3.4 × 1012	4.7 × 1012	9.8 × 1011	1.2 × 1012	6.3 × 1011	8.3 × 1011
x2	Shear stiffness per unit area z2 (N/m3)	1.2 × 1012	1.7 × 1012	3.5 × 1011	4.3 × 1011	2.3 × 1011	3.0 × 1011
x3	Critical normal stress z3 (MPa)	47.2	53.7	21.3	26.8	20.1	22.5
x4	Critical shear stress z4 (MPa)	36.1	39.9	22.2	26.5	17.6	20.2
x5	Bonded disk radius z5 (mm)	0.3	0.8	0.375	1.0	0.45	1.2

Table 3

Factors and levels of the PB test for the grain stalk bonding parameters.

Code	Simulation Parameters	Low Level (−1)	High Level (+1)
x1	Normal stiffness per unit area z1 (N/m3)	6.0 × 109	1.0 × 1010
x2	Shear stiffness per unit area z2 (N/m3)	2.1 × 109	3.6 × 109
x3	Critical normal stress z3 (MPa)	3.5	6.1
x4	Critical shear stress z4 (MPa)	2.9	5.3
x5	Bonded disk radius z5 (mm)	1.875	5

Table 4

PB test plan and results for the stalk bending resistance force.

Order	Normal Stiffness per Unit Area (N/m3)	Shear Stiffness per Unit Area (N/m3)	Critical Normal Stress (MPa)	Critical Shear Stress (MPa)	Bonded Disk Radius (mm)	Bending Resistance Force (N)
						y1	y2	y3
1	1	1	−1	1	1	1.196	2.696	4.860
2	−1	1	1	−1	1	1.202	2.228	5.114
3	1	−1	1	1	−1	0.982	2.540	4.500
4	−1	1	−1	1	1	1.224	2.154	5.066
5	−1	−1	1	−1	1	0.850	2.320	5.062
6	−1	−1	−1	1	−1	0.910	2.244	4.310
7	1	−1	−1	−1	1	0.916	2.410	4.260
8	1	1	−1	−1	−1	0.894	2.140	4.608
9	1	1	1	−1	−1	1.002	2.258	5.004
10	−1	1	1	1	−1	1.158	2.256	5.158
11	1	−1	1	1	1	1.152	2.720	4.804
12	−1	−1	−1	−1	−1	0.772	2.038	4.362

Table 5

Significance analysis of the PB test results for the stalk bending resistance force.

Item	Source	Sum of Squares	Freedom	Mean Square	F-Value	p-Value
y1	Modle	6.32	5	1.26	14.94	0.0025 **
	x1	0.0014	1	0.0014	0.0166	0.9016
	x2	2.49	1	2.49	29.47	0.0016 **
	x3	0.3924	1	0.3924	4.64	0.0747
	x4	2.03	1	2.03	23.94	0.0027 **
	x5	1.41	1	1.41	16.64	0.0065 **
	Residual	0.5077	6	0.0846	/	/
	Cor Total	6.83	11	/	/	/
y2	Modle	11.69	5	2.34	11.25	0.0053 **
	x1	4.84	1	4.84	23.29	0.0029 **
	x2	0.6075	1	0.6075	2.92	0.1381
	x3	0.8533	1	0.8533	4.11	0.0891
	x4	3.08	1	3.08	14.83	0.0085 **
	x5	2.31	1	2.31	11.10	0.0158 *
	Residual	1.25	6	0.2078	/	/
	Cor Total	12.93	11	/	/	/
y3	Modle	28.54	5	5.71	13.70	0.0031 **
	x1	2.24	1	2.24	5.37	0.0597
	x2	13.15	1	13.15	31.55	0.0014 **
	x3	9.86	1	9.86	23.67	0.0028 **
	x4	0.1728	1	0.1728	0.4147	0.5434
	x5	3.12	1	3.12	7.49	0.0339 *
	Residual	2.5	6	0.4167	/	/
	Cor Total	31.04	11	/	/	/

** indicates highly significant (p < 0.01); * indicates significant (p < 0.05).

Table 6

The steepest ascent test plan and results for the stalk bending resistance force.

Item	Parameter	Order
		1	2	3	4	5
Primary stalk	Shear stiffness per unit area (N/m3)	1.20 × 1012	1.325 × 1012	1.45 × 1012	1.575 × 1012	1.70 × 1012
	Critical shear stress (MPa)	36.10	37.05	38.0	38.95	39.90
	Bonded disk radius (mm)	0.30	0.425	0.55	0.675	0.80
	Bending resistance force (N)	0.776	1.097	1.248	1.302	1.363
Secondary stalk	Normal stiffness per unit area (N/m3)	9.8 × 1011	9.9 × 1011	1.0 × 1012	1.1 × 1012	1.2 × 1012
	Critical shear stress (MPa)	22.20	23.275	24.35	25.425	26.50
	Bonded disk radius (mm)	0.375	0.53125	0.6875	0.84375	1.000
	Bending resistance force (N)	1.989	2.133	2.576	2.691	2.724
Tertiary stalk	Shear stiffness per unit area (N/m3)	2.30 × 1011	2.475 × 1011	2.65 × 1011	2.825 × 1011	3.00 × 1011
	Critical normal stress (MPa)	20.10	20.70	21.30	21.90	22.50
	Bonded disk radius (mm)	0.45	0.6375	0.825	1.1025	1.20
	Bending resistance force (N)	4.276	4.353	4.705	5.017	5.189

Table 7

Factors and codes of Box–Behnken test for the stalk bending resistance force.

Item	Parameter	Code
		−1	0	1
primary stalk	Shear stiffness per unit area (N/m3)	1.2 × 1012	1.325 × 1012	1.45 × 1012
	Critical shear stress (MPa)	36.10	37.05	38.00
	Bonded disk radius (mm)	0.3	0.425	0.55
Secondary stalk	Normal stiffness per unit area (N/m3)	9.9 × 1011	1.0 × 1012	1.1 × 1012
	Critical shear stress (MPa)	23.275	24.350	25.425
	Bonded disk radius (mm)	0.53125	0.68750	0.84375
Tertiary stalk	Shear stiffness per unit area (N/m3)	2.475 × 1011	2.65 × 1011	2.825 × 1011
	Critical normal stress (MPa)	20.7	21.3	21.9
	Bonded disk radius (mm)	0.6375	0.8250	1.1025

Table 8

Box–Behnken test plan and results of the stalk bending resistance force.

Test Number	x2/x1/x2	x4/x4/x3	x5/x5/x5	Bending Resistance Force (N)
				y1	y2	y3
1	−1	−1	0	0.796	2.252	4.534
2	1	−1	0	1.068	2.508	4.734
3	−1	1	0	0.972	2.376	4.620
4	1	1	0	1.200	2.664	5.106
5	−1	0	−1	0.692	2.174	4.538
6	1	0	−1	0.922	2.440	4.930
7	−1	0	1	0.780	2.326	4.326
8	1	0	1	1.036	2.576	4.742
9	0	−1	−1	0.786	2.324	4.576
10	0	1	−1	0.938	2.444	4.994
11	0	−1	1	0.794	2.348	4.446
12	0	1	1	1.130	2.702	4.664
13	0	0	0	1.076	2.562	4.910
14	0	0	0	1.094	2.548	4.864
15	0	0	0	1.018	2.574	4.886
16	0	0	0	1.050	2.524	4.782
17	0	0	0	1.006	2.626	4.870

Table 9

ANOVA of bending resistance force for the stalk.

Item	Source	Sum of Squares	Freedom	Mean Square	F-Value	p-Value
y1	Model	0.3302	9	0.0367	26.64	0.0001 **
	x2	0.1215	1	0.1215	88.25	<0.0001 **
	x4	0.0792	1	0.0792	57.51	0.0001 **
	x5	0.0202	1	0.0202	14.67	0.0065 **
	x2x4	0.0005	1	0.0005	0.3515	0.5719
	x2x5	0.0002	1	0.0002	0.1227	0.7364
	x4x5	0.0085	1	0.0085	6.15	0.0423 *
	x22	0.0094	1	0.0094	6.8	0.0351 *
	x42	0.0002	1	0.0002	0.1652	0.6966
	x52	0.0875	1	0.0875	63.53	<0.0001 **
	Residual	0.0096	7	0.0014	/	/
	Lack of fit terms	0.0041	3	0.0014	0.9764	0.4872
	Error	0.0056	4	0.0014	/	/
	Total	0.3399	16	/	/	/
y2	Model	0.3486	9	0.0387	25.65	0.0002 **
	x1	0.1405	1	0.1405	93.0	<0.0001 **
	x4	0.0711	1	0.0711	47.05	0.0002 **
	x5	0.0406	1	0.0406	26.89	0.0013 **
	x1x4	0.0003	1	0.0003	0.1695	0.6929
	x1x5	0.0001	1	0.0001	0.0424	0.8428
	x4x5	0.0137	1	0.0137	9.06	0.0196 *
	x12	0.0389	1	0.0389	25.77	0.0014 **
	x42	0.0018	1	0.0018	1.19	0.3117
	x52	0.0354	1	0.0354	23.42	0.0019 **
	Residual	0.0106	7	0.0015	/	/
	Lack of fit terms	0.0048	3	0.0016	1.11	0.4426
	Error	0.0058	4	0.0014	/	/
	Total	0.3592	16	/	/	/
y3	Model	0.6912	9	0.0768	34.42	<0.0001 **
	x2	0.279	1	0.279	125.05	<0.0001 **
	x3	0.1496	1	0.1496	67.05	<0.0001 **
	x5	0.0925	1	0.0925	41.44	0.0004 **
	x2x3	0.0204	1	0.0204	9.17	0.0192 *
	x2x5	0.0001	1	0.0001	0.0645	0.8068
	x3x5	0.01	1	0.01	4.48	0.072
	x22	0.0237	1	0.0237	10.6	0.0139 *
	x32	0.0064	1	0.0064	2.86	0.1345
	x52	0.0991	1	0.0991	44.44	0.0003 **
	Residual	0.0156	7	0.0022	/	/
	Lack of fit terms	0.0063	3	0.0021	0.8945	0.5172
	Error	0.0093	4	0.0023	/	/
	Total	0.7068	16	/	/	/

** indicates highly significant (p < 0.01); * indicates significant (p < 0.05).

Table 10

PB test plan and results for model threshing rate.

Order	Normal Stiffness per Unit Area (N/m3)	Shear Stiffness per Unit Area (N/m3)	Critical Normal Stress (MPa)	Critical Shear Stress (MPa)	Bonded Disk Radius (mm)	Threshing Rate (%)
1	1	1	−1	1	1	36.96
2	−1	1	1	−1	1	37.58
3	1	−1	1	1	−1	39.43
4	−1	1	−1	1	1	37.16
5	−1	−1	1	−1	1	38.02
6	−1	−1	−1	1	−1	44.92
7	1	−1	−1	−1	1	40.26
8	1	1	−1	−1	−1	39.12
9	1	1	1	−1	−1	35.44
10	−1	1	1	1	−1	36.54
11	1	−1	1	1	1	34.66
12	−1	−1	−1	−1	−1	45.21

Table 11

Significance analysis of the PB test results for model threshing rate.

Source	Sum of Squares	Freedom	Mean Square	F-Value	p-Value
Modle	112.2	5	22.44	12.61	0.0039 **
x1	15.32	1	15.32	8.61	0.0261 *
x2	32.34	1	32.34	18.17	0.0053 **
x3	40.19	1	40.19	22.58	0.0032 **
x4	2.96	1	2.96	1.66	0.2446
x5	21.39	1	21.39	12.02	0.0134 *
Residual	10.68	6	1.78	/	/
Cor Total	122.87	11	/	/	/

** indicates highly significant (p < 0.01); * indicates significant (p < 0.05).

Table 12

The steepest ascent test plan and the results for the model threshing rate.

Order	Normal Stiffness per Unit Area (N/m3)	Shear Stiffness per Unit Area (N/m3)	Critical Normal Stress (MPa)	Bonded Disk Radius (mm)	Threshing Rate (%)
1	1.0 × 1010	3.6 × 109	6.1	5	35.33
2	9.0 × 109	3.225 × 109	5.45	4.21875	36.46
3	8.0 × 109	2.85 × 109	4.8	3.4375	38.14
4	7.0 × 109	2.475 × 109	4.15	2.65625	41.22
5	6.0 × 109	2.1 × 109	3.5	1.875	43.18

Table 13

Factors and codes of Box–Behnken test for model threshing rate.

Code	Normal Stiffness per Unit Area (N/m3)	Shear Stiffness per Unit Area (N/m3)	Critical Normal Stress (MPa)	Bonded Disk Radius (mm)
−1	7.0 × 109	2.475 × 109	4.15	2.625
0	8.0 × 109	2.85 × 109	4.80	2.75
+1	9.0 × 109	3.225 × 109	5.45	2.875

Table 14

Box–Behnken test plan and results for model threshing rate.

Test Number	x1	x2	x3	x5	y4 (%)
1	−1	−1	0	0	40.49
2	1	−1	0	0	38.14
3	−1	1	0	0	35.85
4	1	1	0	0	36.72
5	0	0	−1	−1	38.42
6	0	0	1	−1	35.83
7	0	0	−1	1	36.22
8	0	0	1	1	34.05
9	−1	0	0	−1	38.76
10	1	0	0	−1	37.43
11	−1	0	0	1	37.37
12	1	0	0	1	34.93
13	0	−1	−1	0	39.21
14	0	1	−1	0	35.74
15	0	−1	1	0	35.32
16	0	1	1	0	34.81
17	−1	0	−1	0	40.38
18	1	0	−1	0	37.57
19	−1	0	1	0	36.11
20	1	0	1	0	36.04
21	0	−1	0	−1	37.33
22	0	1	0	−1	36.74
23	0	−1	0	1	36.64
24	0	1	0	1	33.45
25	0	0	0	0	37.71
26	0	0	0	0	37.25
27	0	0	0	0	36.96
28	0	0	0	0	36.94
29	0	0	0	0	36.43
30	0	0	0	0	36.34
31	0	0	0	0	37.44

Table 15

ANOVA of bending resistance force for model threshing rate.

Source	Sum of Squares	Freedom	Mean Square	F-Value	p-Value
Model	73.19	14	5.23	24.16	<0.0001 **
x1	5.51	1	5.51	25.46	0.0001 **
x2	15.92	1	15.92	73.56	<0.0001 **
x3	19.71	1	19.71	91.11	<0.0001 **
x5	11.7	1	11.7	54.09	<0.0001 **
x1x2	2.59	1	2.59	11.98	0.0032 **
x1x3	1.88	1	1.88	8.68	0.0095 **
x1x4	0.308	1	0.308	1.42	0.2502
x2x3	2.19	1	2.19	10.12	0.0058 **
x2x4	1.69	1	1.69	7.81	0.013 *
x3x4	0.0441	1	0.0441	0.2038	0.6577
x12	5.84	1	5.84	27.01	<0.0001 **
x22	0.5005	1	0.5005	2.31	0.1478
x32	0.9116	1	0.9116	4.21	0.0568
x52	3.24	1	3.24	14.98	0.0014 **
Residual	3.46	16	0.2164	/	/
Lack of fit terms	1.94	10	0.1936	0.7618	0.665
Error	1.53	6	0.2542	/	/
Total	76.65	30	/	/	/

** indicates highly significant (p < 0.01); * indicates significant (p < 0.05).

Appendix A

As illustrated in Figure A1a, a rectangular coordinate system O-XYZ was established by taking the center of the bottom circle of the third internode as the coordinate origin O. From this origin O, the particles of the third internode model of rice were stacked upward along the z-axis. The coordinates of these particles can be expressed using Equation (A1).(A1)$$\left(X_{s_3}^k{,Y_{s_3}^k,Z_{s_3}^k)=(\mathrm{0,0},2(k-1)R}_{S_3}\right)$$ where X_S3^k, Y_S3^k, and Z_S3^k are the three-axis coordinates of the kth particle in the tertiary stalk, k = 1…384, and R_S3 is the radius of the tertiary stalk particle.

The coordinates of the secondary stalk particles are superimposed on the coordinates of the tertiary stalk particles.(A2)$$\left(X_{s_2}^k{,Y_{s_2}^k,Z_{s_2}^k)=(\mathrm{0,0},2N_{S_3}^{k-1}{R}_{S_3}+2(k-1)R}_{S_2}\right)$$ where X_S2^k, Y_S2^k, and Z_S2^k are the three-axis coordinates of the kth particle in the secondary stalk, k = 1…502, and R_S2 is the radius of the secondary stalk particle.

The coordinates of the different stalk particles were accumulated on the previous cumulative stalk particle model, so the coordinates of the primary stalk can be expressed as:(A3)$$\left(X_{s_1}^k,Y_{s_1}^k,Z_{s_1}^k)=(\mathrm{0,0},2N_{S_3}^{j-1}{R}_{S_3}+2N_{S_2}^{j-1}{R}_{S_2}+2(k-1)R_{S_1}\right)$$ where X_S1^k, Y_S1^k, and Z_S1^k are the three-axis coordinates of the kth particle in the primary stalk, k = 1…607, R_S1 is the radius of the primary stalk particle, N_S2^j and N_S3^j represents the number of particles in the secondary stalk and the tertiary stalk.

Simplified to Equation (A4):(A4)$$\left(X_{s_1}^k,Y_{s_1}^k,Z_{s_1}^k)=(\mathrm{0,0},2\sum _{i=2}^3{N}_{S_n}^{j-1}{R}_{S_n}+2(k-1)R_{S_1}\right)$$

Similarly, the coordinates of the cob particles can be expressed as Equation (A5):(A5)$$\left(X_r^k,Y_r^k,Z_r^k)=(\mathrm{0,0},2\sum _{i=1}^n{N}_{S_n}^{j-1}{R}_{S_n}+2(k-1)R_r\right)$$ where X_r^k, Y_r^k, and Z_r^k are the three-axis coordinate of the kth particle in the rachis, k = 1…64, R_r is the radius of the rachis particle, and n is the number of internodes, n = 1,2,3. N_Sn^j indicates the number of particles in the nth stalk, j = 607,502,384, and R_Sn indicates the radius of the nth stalk particle.

Taking the Y-Z plane as an example, the side primary branches were symmetrically distributed around the z axis, and the included angle is the growth angle ±α. Taking the center of each particle, where the primary branch contacts the rachis as the coordinate origin O₁, the z axis is the starting boundary, and rotating around the x axis by ±α, the local coordinate system Y₁O₁Z₁ is established, as shown in Figure A1b. The coordinate position of the primary branch particle is determined by homogeneous coordinate transformation.

The global coordinates of the particle where the primary branch contacts the rachis can be expressed as Equation (A6):(A6)$$\left(X_p^k,Y_p^k,Z_p^k\right)=\left(0,{\left(-1\right)}^{k-1}\left(R_r+R_P\right),2\sum _{i=1}^n{N}_{S_n}^{j-1}{R}_{S_n}+\left[{\displaystyle \frac{k-1}{2}}\right]L_1+L\right)$$ where X_p^k, Y_p^k, and Z_p^k are the global coordinates of the kth particle in contact with the rachis, k = 1…6. R_p is the radius of the primary branch particle, L₁ is the distance between adjacent primary branches along the z-axis, L is the distance between the first particle and the top particle of the primary stalk, and $\left[{\displaystyle \frac{k-1}{2}}\right]$ is the rounding function.

The coordinates of the particles of the primary branch in the local coordinate system Y₁O₁Z₁ can be expressed as Equation (A7):(A7)$$\left(x_p^k{,y_p^k,z_p^k)=(\mathrm{0,0},2(k-1)R}_p\right)$$ where x_p^k, y_p^k, and z_p^k are the three-axis coordinates of the kth particle in the local coordinate system Y₁O₁Z₁ where the primary branch and the rachis come into contact. k = 1…n_p, R_p is the radius of the primary branch particle.

The rotary and translational matrix C₁ between the global coordinate system and the local coordinate system Y₁O₁Z₁ can be expressed as Equation (A8):(A8)$$C_1=R_1{T}_1=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& cos\alpha & {\left(-1\right)}^{k}sin\alpha & 0\\ 0& {\left(-1\right)}^{k-1}sin\alpha & cos\alpha & 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& {\left(-1\right)}^{k-1}\left(R_r+R_P\right)& 2{\displaystyle \sum _{i=1}^n}{N}_{S_n}^{j-1}{R}_{S_n}+L& 1\end{array}\right]$$ where R₁ is the rotation transformation matrix under the rotation angle of ±α and T₁ is the translation transformation matrix of local coordinate system 1.

From this, the global coordinates of any side primary branch particle in the Y-Z plane can be obtained.(A9)$$\left(X_p^k,Y_p^k,Z_p^k\right)=R_1{T}_1\left(x_p^k,y_p^k,z_p^k\right)$$

Similarly, we can obtain the grain stalk coordinates under the local coordinates Y₁O₁Z₁.(A10)$$(x_{GS}^k{y}_{GS}^k,z_{GS}^k)=(0,{\left(-1\right)}^{k-1}{(R}_P+R_{GS}),\left[{\displaystyle \frac{k-1}{2}}\right]L_2+L^{\prime })$$ where x_GS^k, y_GS^k, and z_GS^k are the coordinates of the kth grain stalk particle in the local coordinate system 1 and k = 1…n_GS. R_GS is the radius of the grain stalk particle, L′ is the distance between the first grain stalk and the starting end of the primary branch, and L₂ is the distance between adjacent grain stalk particles.

After rotary and translational coordinate transformation, the global coordinates of the peduncle on the side primary branch can be expressed as Equation (A11):(A11)$$\left(X_{gs}^k,Y_{gs}^k,Z_{gs}^k\right)=C_1\left(x_{gs}^k{y}_{gs}^k,z_{gs}^k\right)$$

Rice grains are connected by a primary branch of the grain stalk, and the included angle is the growth angle ±β. With the center of the grain stalk and the grain as the coordinate origin O₂ and the growth direction of the lateral primary branch as the Z₂ direction, a local coordinate system Y₂O₂Z₂ is established, as shown in Figure A1c. With the center of the grain-end grain as the local coordinate origin O₃, the Z₂ axis as the starting boundary, and the rotation around the x-axis ±β, a local coordinate system Y₃O₃Z₃ is established, as shown in Figure A1d. The coordinate position of the grain is determined by rotary and translational coordinate transformation.

The particle coordinates of the grain in the local coordinate system Y₂O₂Z₂ can be expressed as Equation (A12):(A12)$$\left(x_G^k{y}_G^k,z_G^k\right)=(0,{\left(-1\right)}^{k-1}{(R}_G+R_{GS}),\left[{\displaystyle \frac{k-1}{2}}\right]L_3)$$ where x_G^k, y_G^k, and z_G^k are the coordinates of the kth handle in the local coordinate system Y₂O₂Z₂, k = 1…n_G. R_G is the radius of the top of the grain. L₃ is the distance between adjacent grains.

The rotary and translational transformation matrix between the local coordinate system Y₃O₃Z₃ and the local coordinate system Y₂O₂Z₂ can be expressed as Equation (A13):(A13)$$C_2=R_2{T}_2$$(A14)$$T_2=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& {\left(-1\right)}^{k-1}{(R}_P+R_{GS})& L^{\prime }& 1\end{array}\right]$$(A15)$$R_2=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& cos\beta & {\left(-1\right)}^{k}sin\beta & 0\\ 0& {\left(-1\right)}^{k-1}sin\beta & cos\beta & 0\\ 0& 0& 0& 1\end{array}\right]$$ where R₂ is the rotation transformation matrix under the rotation angle of ±β and T₂ is the translation transformation matrix of the local coordinate system Y₃O₃Z₃.

The rotary and translational transformation matrix between the global coordinate system and the local coordinate system Y₃O₃Z₃ can be expressed as Equation (A16):(A16)$$C_3=R_2{R}_1{T}_2{T}_1$$

This allows us to obtain the global coordinates of the grain on any side primary branch in the Y-Z plane.(A17)$$\left(X_G^k,Y_G^k,Z_G^k\right)=C_3\left(x_G^k{y}_G^k,z_G^k\right)$$

The coordinates of the top primary branch particles in the global coordinate system can be expressed as Equation (A18):(A18)$$\left(X_{tp}^k,Y_{tp}^k,Z_{tp}^k\right)=\left(0,{\left(-1\right)}^{k-1}\left(R_r+R_P\right),2N_r^{k-1}{R}_r+2\sum _{i=1}^n{N}_{S_n}^{j-1}{R}_{S_n}+2\left(k-1\right)R_p\right)$$

The coordinates of the grain on the top primary branch in the global coordinate system can be expressed as Equation (A19):(A19)$$\left(X_{tgs}^k,Y_{tgs}^k,Z_{tgs}^k\right)=\left(0,{\left(-1\right)}^{k-1}\left(R_r+{2R}_P+R_{gs}\right),2N_r^{k-1}{R}_r+2\sum _{i=1}^n{N}_{S_n}^{j-1}{R}_{S_n}+\left[{\displaystyle \frac{k-1}{2}}\right]L_2+L^{\prime }\right)$$

Grains on the top primary branch are symmetrically distributed around the z axis, and the included angle is the growth angle ±β. Take the center of the grain on the top primary branch as the coordinate origin, the z axis as the starting boundary, rotate around the x axis ±β, and establish the local coordinate system Y₄O₄Z₄, as shown in Figure A1e. The grain coordinates of the grain in the local coordinate system Y₄O₄Z₄ can be expressed as Equation (A20):(A20)$$\left(x_{TG}^k{y}_{TG}^k,z_{TG}^k\right)=(0,{\left(-1\right)}^{k-1}{(R}_G+R_{gs}),\left[{\displaystyle \frac{k-1}{2}}\right]L_3)$$ where x_TG^k, y_TG^k, and z_TG^k are the coordinates of the kth grain stalk in the local coordinate system Y₄O₄Z₄ and k = 1…n_TG. R_TG is the radius of the top of the grain. L₃ is the distance between adjacent grains.

The rotary and translational transformation matrix between the global coordinate system and the local coordinate system Y₄O₄Z₄ can be expressed as Equation (A21):(A21)$$C_4=R_2{T}_3=\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& cos\beta & {\left(-1\right)}^{k}sin\beta & 0\\ 0& {\left(-1\right)}^{k-1}sin\beta & cos\beta & 0\\ 0& 0& 0& 1\end{array}\right]\left[\begin{array}{cccc}1& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& {\left(-1\right)}^{k-1}\left(R_r+{2R}_P+R_{gs}\right)& 2N_r^{j-1}{R}_r+2{\displaystyle \sum _{i=1}^n}{N}_{S_n}^{j-1}{R}_{S_n}+L^{\prime }& 1\end{array}\right]$$ where T₃ is the translation transformation matrix of the local coordinate system Y₄O₄Z₄.

From this, we can obtain the global coordinates of the grain on the top primary branch.(A22)$$\left(X_{TG}^k,Y_{TG}^k,Z_{TG}^k\right)=C_4\left(x_{TG}^k{y}_{TG}^k,z_{TG}^k\right)$$

Figure A1. Coordinate system of the rice discrete element model. (a) Original coordinate system of the discrete element model O-XYZ. (b) Rotary and translational coordinate system of the lateral primary branch Y1O1Z1. (c) Translational coordinate system of the grain stalk Y2O2Z2. (d) Rotary coordinate system of the grain Y3O3Z3. (e) Translational coordinate system of the top primary branch Y4O4Z4.

Appendix B

Table A1

Rice biological characteristic parameters.

Plant Characteristics		Minimum Value	Maximum Value	Mean Value	Standard Deviation
Weight	Thousand-grain weight (g)	30.97	32.85	31.85	0.652
Moisture content	Grain moisture content (%)	26.8	29.3	27.94	0.823
	Stalk moisture content (%)	51.4	55.9	53.60	1.487
	Stalk moisture content (%)	68.5	73.3	71.18	1.675
Density	Grain (kg/m3)	853	1149	1004.95	99.755
	Branch (kg/m3)	165	243	209.8	25.289
	Stalk (kg/m3)	186	268	224.60	28.196
Length	Grain stalk length (cm)	0.1	0.5	0.22	0.114
	Primary branch length (cm)	4.2	7.6	6.43	0.966
	Rachis length (cm)	6.2	9.6	8.46	1.176
	Primary stalk length (cm)	17.3	29.0	24.29	3.005
	Secondary stalk length (cm)	16.5	29.9	25.12	3.141
	Tertiary stalk length (cm)	17.7	29.5	23.01	3.602
Diameter	Diameter of grain stalk (mm)	0.16	0.97	0.59	0.280
	Diameter of primary branch (mm)	0.67	0.96	0.82	0.093
	Diameter of rachis (mm)	1.46	2.59	1.99	0.349
	Outer diameter of primary stalk (mm)	1.71	2.22	1.96	0.158
	Outer diameter of secondary stalk (mm)	2.80	3.44	3.08	0.202
	Outer diameter of tertiary stalk (mm)	3.63	4.82	4.02	0.256
	Inner diameter of primary stalk (mm)	1.03	1.49	1.28	0.146
	Inner diameter of secondary stalk (mm)	1.66	2.49	2.03	0.229
	Inner diameter of tertiary stalk (mm)	2.45	3.70	2.77	0.213
Grain size	Grain length (mm)	7.96	8.66	8.29	0.203
	Grain width (mm)	3.49	4.05	3.89	0.137
	Grain thickness (mm)	2.55	2.75	2.65	0.058

Table A2

Methods for measuring physical parameters of rice.

Measurement Items	Determination Method	Calculation
Moisture content	The moisture content of various parts of a rice plant exhibits considerable variability, necessitating the measurement of moisture levels in different plant components. The rachis, primary branches, and secondary branches were collectively referred to as branches. A drying method was employed to assess the moisture content of rice grains, branches, and stalks. The testing equipment included an electric heating air-drying oven (model: 101-00S, accuracy: ±0.1 °C, Zhejiang Tianyu intelligent Technology Co., Ltd., Wenzhou, China) and an electronic analytical balance (model: YH-M, range 600 g, accuracy: ±0.01 g, Yingheng Electronic Technology Co., Ltd, Huizhou, China). Initially, the weight of the stalk was measured using the electronic scale and recorded. The stalk was then placed in a blast drying oven at 150 °C for 12 h, or until its mass stabilized. Afterward, the dried mass was weighed and recorded. The moisture content of the test sample was calculated according to Equation (1). This process was repeated 10 times for each group to determine the average value and standard deviation. The same methodology was applied to measure the moisture content of the rice grains and branches.	$\begin{array}{c}w={\displaystyle \frac{m-m^{\prime }}{m}}\times 100\%\end{array}$where m represents the initial mass of the sample, kg, and m′ represents the mass of the sample after drying, kg
Density	Density was estimated using the immersion substitution method. Given the complexity of the rice structure, primary branches, secondary branches, and the rachis were collectively classified under the stalk category, with the densities of grains and stalks measured separately. A sample was weighed and immersed in a measuring cylinder with a volume of 100 mL, which initially contained 50 mL of water. The density of the sample was then calculated using Equation (2). This process was repeated 10 times, and the average density was recorded.	$\rho ={\displaystyle \frac{m}{V_2-V_1}}$where m represents the mass of the sample, kg; V1 represents the volume in the measuring cylinder before adding the sample, m3; and V2 represents the volume in the measuring cylinder before adding the sample, m3.
Coefficient of static friction	The static friction coefficient was measured using an angle adjustment platform (model: MGO-000-0032, adjustment range 0–90°, GWOOD Inc., Jinzhou, China). Taking grain–steel as an example, the steel plate was installed on the adjustment platform, which was positioned parallel to the ground. Grains were then placed on the angle adjustment platform. The angle adjustment knob was rotated to gradually increase the angle of the platform. When the grains began to move or exhibited a tendency to move, the rotation was halted, and the angle sensor was read. This process was repeated ten times to calculate the average and standard deviation. The same procedure was employed to measure other coefficients of static friction.	$\begin{array}{c}{\mu }_1=\mathit{tan}\lambda \end{array}$where μ1 is the static friction coefficient and λ is the angle at which the motion trend is generated.
Coefficient of dynamic friction	The dynamic friction coefficient measurement test was conducted using a high-speed camera (model: FASTCAM-MH6, minimum resolution 1920 × 1400, Photron Inc., Tokyo, Japan). Taking grain–steel as an example, a steel plate was installed on the angle adjustment platform, and the angle was fixed at 45 degrees. At a height of 40 mm, the grains rolled down the slope to the bottom under the influence of gravity. The high-speed camera recorded the motion of the grains as they traveled from the top to the bottom. During this process, the average speed of the grains was calculated by analyzing five frames captured before and after reaching the bottom, utilizing the frame rate of the camera and grid paper coordinates; this value was recorded as the instantaneous speed at the bottom. This procedure was repeated 10 times to obtain an average value. The same methodology was employed to measure other dynamic friction coefficients.	$\begin{array}{c}g{h}_1={\mu }_{2}g{h}_1\mathit{cot}\varphi +{\displaystyle \frac{1}{2}}{\nu }^2\end{array}$where μ2 is the coefficient of kinetic friction; h1 is the height of the slope, m; η is the angle of the slope; v is the speed at which the object rolls to the bottom, m/s; and g is the acceleration due to gravity, m/s2.
Coefficient of restitution	Using grain–steel as an example, a steel plate was installed on the angle adjustment platform. The angle adjustment platform was then adjusted to ensure it was parallel to the ground. A free-fall platform was positioned 100 mm above the surface of the steel plate, onto which grains were placed to fall. Upon reaching the bottom of the steel plate, the grains would rebound. A high-speed camera was employed to record the maximum height achieved by the grains after their contact with the steel plate. This process was repeated ten times to obtain an average value. Similarly, the recovery coefficients for other materials can be measured.	$\begin{array}{c}e=\sqrt{{\displaystyle \frac{h_3}{h_2}}}\end{array}$where e is the coefficient of restitution; h2 is the free fall height, m; and h3 is the maximum rebound height, m.

Table A3

Contact parameters.

Parameters		Minimum Value	Maximum Value	Mean Value	Standard Deviation	Source
Coefficient of static friction	Grain–Grain	0.719	1.510	0.898	0.121	Grade test
	Stalk–Stalk	0.376	0.613	0.519	0.079	Grade test
	Stalk–Grain	0.449	0.695	0.582	0.069	Grade test
	Steel–Grain	0.529	0.713	0.654	0.057	Grade test
	Steel–Stalk	0.414	0.603	0.488	0.058	Grade test
	Grain–Acrylic	0.392	0.579	0.497	0.057	Grade test
	Stalk–Acrylic	0.316	0.548	0.451	0.078	Grade test
Coefficient of rolling friction	Grain–Grain	0.0256	0.0681	0.050	0.014	Slope rolling test
	Stalk–Stalk	0.0119	0.0428	0.021	0.002	Slope rolling test
	Stalk–Grain	0.0214	0.0331	0.027	0.004	Slope rolling test
	Steel–Grain	0.0139	0.0366	0.028	0.006	Slope rolling test
	Steel–Stalk	0.0199	0.0343	0.025	0.004	Slope rolling test
	Grain–Acrylic	0.0184	0.0317	0.025	0.004	Slope rolling test
	Stalk–Acrylic	0.0108	0.0205	0.016	0.003	Slope rolling test
Coefficient of restitution	Grain–Grain	0.316	0.507	0.383	0.064	Drop test
	Stalk–Stalk	0.458	0.806	0.615	0.107	Drop test
	Stalk–Grain	0.346	0.640	0.503	0.090	Drop test
	Steel–Grain	0.361	0.794	0.5453	0.131	Drop test
	Steel–Stalk	0.436	0.663	0.552	0.063	Drop test
	Grain–Acrylic	0.433	0.572	0.511	0.043	Drop test
	Stalk–Acrylic	0.517	0.694	0.615	0.056	Drop test

Table A4

Normal stiffness of different structures of rice plant.

Parameters	Maximum Tensile Force (N/mm)
	Minimum Value	Maximum Value	Mean Value	Standard Deviation
Grain stalk	1.18	1.96	1.67	0.232
Primary branches	2.24	4.11	3.26	0.613
Rachis	19.92	24.52	22.55	1.521
Primary stalks	20.68	28.56	25.31	2.324
Secondary stalks	26.72	31.90	29.72	1.688
Tertiary stalks	33.69	44.61	39.43	3.509
Primary branches–rachis	1.65	2.36	2.12	0.230
Rachis–primary stalk	22.93	28.90	26.39	1.900
Primary stalks–secondary stalks	37.43	46.02	42.34	3.026
Secondary stalks–tertiary stalks	45.33	55.79	50.91	3.172

Table A5

Maximum tensile forces of different structures of rice plant.

Parameters	Maximum Tensile Force (N)
	Minimum Value	Maximum Value	Mean Value	Standard Deviation
Grain stalk	0.95	1.65	1.29	0.227
Primary branches	4.50	6.48	5.44	0.655
Rachis	31.26	38.06	35.17	2.433
Primary stalks	50.06	56.95	52.89	2.0587
Secondary stalks	60.21	75.76	68.59	5.347
Tertiary stalks	79.59	89.10	84.66	3.195
Primary branches–rachis	4.50	6.12	5.58	0.398
Rachis–primary stalk	38.99	48.90	42.69	2.8438
Primary stalks–secondary stalks	96.34	109.04	102.67	4.398
Secondary stalks–tertiary stalks	101.86	124.18	113.54	7.418

Table A6

Maximum shear forces of different structures of rice plant.

Parameters	Maximum Tensile Force (N)
	Minimum Value	Maximum Value	Mean Value	Standard Deviation
Grain stalk	0.78	1.43	1.12	0.2108
Primary branches	21.55	26.46	23.78	1.599
Rachis	27.54	34.66	31.51	2.317
Primary stalks	38.28	42.31	39.98	1.280
Secondary stalks	62.75	74.91	68.80	3.968
Tertiary stalks	69.69	79.99	75.08	3.278
Primary branches–rachis	3.77	5.39	4.61	0.507
Rachis–primary stalk	42.09	60.04	50.36	5.716
Primary stalks–secondary stalks	79.40	106.92	95.22	8.811
Secondary stalks–tertiary stalks	114.41	179.99	150.32	22.063

Table A7

Bond parameters of rice plants in different structures.

Parameter	Normal Stiffness per Unit Area (N/m3)	Shear Stiffness per Unit Area (N/m3)	Critical Normal Stress (MPa)	Critical Shear Stress (MPa)
Grain stalk	6.0 × 109–1.0 × 1010	2.1 × 109–3.6 × 109	3.5–6.1	2.9–5.3
Primary branches	1.8 × 1011–3.3 × 1011	6.4 × 1010–1.2 × 1011	8.6–12.4	41.2–50.6
Rachis	2.6 × 1011–3.2 × 1011	9.3 × 1010–1.1 × 1011	10.1–12.3	8.9–11.2
Primary stalks	3.4 × 1012–4.7 × 1012	1.2 × 1012–1.7 × 1012	47.2–53.7	36.1–39.9
Secondary stalks	9.8 × 1011–1.2 × 1012	3.5 × 1011–4.3 × 1011	21.3–26.8	22.2–26.5
Tertiary stalks	6.3 × 1011–8.3 × 1011	2.3 × 1011–3.0 × 1011	20.1–22.5	17.6–20.2
Primary branches–rachis	8.4 × 109–1.2 × 1010	3.0 × 109–4.3 × 109	8.6–11.7	7.2–10.3
Rachis–primary stalk	1.72 × 1012–2.17 × 1012	6.1 × 1011–7.8 × 1011	12.6–15.8	13.6–19.4
Primary stalks–secondary stalks	6.6 × 1011–8.1 × 1011	2.4 × 1011–2.9 × 1011	9.1–10.3	7.5–10.1
Secondary stalks–tertiary stalks	5.4 × 1011–6.6 × 1011	1.9 × 1011–2.0 × 1011	7.3–8.9	8.2–12.9

Table A8

Threshing rate at different impact angles.

Impact Angle	Minimum Value	Maximum Value	Mean Value	Standard Deviation
15°	12.13	18.46	14.76	2.16
30°	22.79	31.07	26.37	2.53
45°	37.21	42.45	39.71	1.81
60°	55.16	63.23	59.71	2.31
75°	61.28	69.03	64.81	2.81
90°	72.57	80.26	76.22	1.98

References

1. Zhao, Z.; Huang, H.; Yin, J.; Yang, S. Dynamic analysis and reliability design of round baler feeding device for rice straw harvest. Biosyst. Eng.; 2018; 174, pp. 10-19. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2018.06.014]

2. Zami, M.A.; Hossain, M.A.; Sayed, M.A.; Biswas, B.K.; Hossain, M.A. Performance evaluation of the BRRI reaper and Chinese reaper compared to manual harvesting of rice (Oryza sativa L.). Agriculturists; 2014; 12, pp. 142-150. [DOI: https://dx.doi.org/10.3329/agric.v12i2.21743]

3. Tang, H.; Xu, C.; Zhao, J.; Wang, J. Stripping mechanism and loss characteristics of a stripping-prior-to-cutting header for rice harvesting based on CFD-DEM simulations and bench experiments. Biosyst. Eng.; 2023; 229, pp. 116-136. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2023.03.023]

4. Tang, Z.; Li, Y.; Li, X.; Xu, T. Structural damage modes for rice stalks undergoing threshing. Biosyst. Eng.; 2019; 186, pp. 323-336. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2019.08.005]

5. Sirikun, C.; Samseemoung, G.; Soni, P.; Langkapin, J.; Srinonchat, J. A grain yield sensor for yield mapping with local rice combine harvester. Agriculture; 2021; 11, 897. [DOI: https://dx.doi.org/10.3390/agriculture11090897]

6. Li, Y.; Wang, X.; Xu, L. Threshing injury to rice grain based on the energy conservation. J. Mech. Eng.; 2007; 43, pp. 160-164. [DOI: https://dx.doi.org/10.3901/JME.2007.03.160]

7. Zeng, Z.; Ma, X.; Cao, X.; Li, Z.; Wang, X. Critical review of applications of discrete element method in agricultural engineering. Trans. Chin. Soc. Agric. Mach.; 2021; 52, pp. 1-20. [DOI: https://dx.doi.org/10.6041/j.issn.1000-1298.2021.04.001]

8. Du, C.; Han, D.; Song, Z.; Chen, Y.; Chen, X.; Wang, X. Calibration of contact parameters for complex shaped fruits based on discrete element method: The case of pod pepper (Capsicum annuum). Biosyst. Eng.; 2023; 226, pp. 43-54. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2022.12.005]

9. Wu, J.; Tang, Q.; Mu, S.; Jiang, L.; Hu, Z. Test and optimization of oilseed rape (Brassica napus L.) threshing device based on DEM. Agriculture; 2022; 12, 1580. [DOI: https://dx.doi.org/10.3390/agriculture12101580]

10. Tian, Y.; Zeng, Z.; Gong, H.; Zhou, Y.; Qi, L.; Zhen, W. Simulation of tensile behavior of tobacco leaf using the discrete element method (DEM). Comput. Electron. Agric.; 2023; 205, 107570. [DOI: https://dx.doi.org/10.1016/j.compag.2022.107570]

11. Yu, Q.; Liu, Y.; Chen, X.; Sun, K.; Lai, Q. Calibration and experiment of simulation parameters for Panax notoginseng seeds based on DEM. Trans. CSAM; 2020; 51, pp. 123-132. [DOI: https://dx.doi.org/10.6041/j.issn.1000-1298.2020.02.014]

12. Tang, Z.; Zhang, B.; Wang, B.; Wang, M.; Chen, H.; Li, Y. Breaking paths of rice stalks during threshing. Biosyst. Eng.; 2021; 204, pp. 346-357. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2021.02.008]

13. Lenaerts, B.; Aertsen, T.; Tijskens, E.; De Ketelaere, B.; Ramon, H.; De Baerdemaeker, J.; Saeys, W. Simulation of grain–straw separation by Discrete Element Modeling with bendable straw particles. Comput. Electron. Agric.; 2014; 101, pp. 24-33. [DOI: https://dx.doi.org/10.1016/j.compag.2013.12.002]

14. Leblicq, T.; Vanmaercke, S.; Ramon, H.; Saeys, W. Mechanical analysis of the bending behaviour of plant stems. Biosyst. Eng.; 2015; 129, pp. 87-99. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2014.09.016]

15. Leblicq, T.; Smeets, B.; Vanmaercke, S.; Ramon, H.; Saeys, W. A discrete element approach for modelling bendable crop stems. Comput. Electron. Agric.; 2016; 124, pp. 141-149. [DOI: https://dx.doi.org/10.1016/j.compag.2016.03.022]

16. Mao, H.; Wang, Q.; Li, Q. Modelling and simulation of the straw-grain separation process based on a discrete element model with flexible hollow cylindrical bonds. Comput. Electron. Agric.; 2020; 170, 105229. [DOI: https://dx.doi.org/10.1016/j.compag.2020.105229]

17. Shi, Y.; Jiang, Y.; Wang, X.; Thuy, N.T.D.; Yu, H. A mechanical model of single wheat straw with failure characteristics based on discrete element method. Biosyst. Eng.; 2023; 230, pp. 1-15. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2023.03.017]

18. Schramm, M.; Tekeste, M.Z.; Plouffe, C.; Harby, D. Estimating bond damping and bond Young’s modulus for a flexible wheat straw discrete element method model. Biosyst. Eng.; 2019; 186, pp. 349-355. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2019.08.003]

19. Su, Z.; Li, Y.; Dong, Y.; Tang, Z.; Liang, Z. Simulation of rice threshing performance with concentric and non-concentric threshing gaps. Biosyst. Eng.; 2020; 197, pp. 270-284. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2020.05.020]

20. Liu, Y.; Li, Y.; Zhang, T.; Huang, M. Effect of concentric and non-concentric threshing gaps on damage of rice straw during threshing for combine harvester. Biosyst. Eng.; 2022; 219, pp. 1-10. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2022.04.013]

21. Liu, F.; Zhang, J.; Chen, J. Modeling of flexible wheat straw by discrete element method and its parameter calibration. Int. J. Agric. Biol. Eng.; 2018; 11, pp. 42-46. [DOI: https://dx.doi.org/10.25165/j.ijabe.20181103.3381]

22. Wang, J.; Xu, C.; Xu, W.; Fu, Z.; Wang, Q.; Tang, H. Discrete element method simulation of rice grain motion during discharge with an auger operated at various inclinations. Biosyst. Eng.; 2022; 223, pp. 97-115. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2022.08.020]

23. Tang, H.; Xu, W.; Zhao, J.; Xu, C.; Wang, J. Comparison of rice straw compression characteristics in vibration mode based on discrete element method. Biosyst. Eng.; 2023; 230, pp. 191-204. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2023.04.009]

24. Tang, H.; Zhu, G.; Xu, W.; Xu, C.; Wang, J. Discrete element method simulation of rice grains impact fracture characteristics. Biosyst. Eng.; 2024; 237, pp. 50-70. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2023.11.011]

25. Guo, J.; Karkee, M.; Yang, Z.; Fu, H.; Li, J.; Jiang, Y.; Duan, J. Discrete element modeling and physical experiment research on the biomechanical properties of banana bunch stalk for postharvest machine development. Comput. Electron. Agric.; 2021; 188, 106308. [DOI: https://dx.doi.org/10.1016/j.compag.2021.106308]

26. Sadrmanesh, V.; Chen, Y. Simulation of tensile behavior of plant fibers using the discrete element method (DEM). Compos. Part A Appl. Sci. Manuf.; 2018; 114, pp. 196-203. [DOI: https://dx.doi.org/10.1016/j.compositesa.2018.08.023]

27. Wang, Y.; Zhang, Y.; Yang, Y.; Zhao, H.; Yang, C.; He, Y.; Xu, H. Discrete element modelling of citrus fruit stalks and its verification. Biosyst. Eng.; 2020; 200, pp. 400-414. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2020.10.020]

28. Zhao, W.; Chen, M.; Xie, J.; Cao, S.; Wu, A.; Wang, Z. Discrete element modeling and physical experiment research on the biomechanical properties of cotton stalk. Comput. Electron. Agric.; 2023; 204, 107502. [DOI: https://dx.doi.org/10.1016/j.compag.2022.107502]

29. Ding, X.; Wang, B.; He, Z.; Shi, Y.; Li, K.; Cui, Y.; Yang, Q. Fast and precise DEM parameter calibration for Cucurbita ficifolia seeds. Biosyst. Eng.; 2023; 236, pp. 258-276. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2023.11.004]

30. Dong, J.; Zhang, D.; Yang, L.; Cui, T.; Zhang, K.; He, X.; Jing, M. Discrete element method optimisation of threshing components to reduce maize kernel damage at high moisture content. Biosyst. Eng.; 2023; 233, pp. 221-240. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2023.08.005]

31. Zhang, D.; Yi, S.; Zhang, J.; Bao, Y. Establishment of millet threshing and separating model and optimization of harvester parameters. Alex. Eng. J.; 2022; 61, pp. 11251-11265. [DOI: https://dx.doi.org/10.1016/j.aej.2022.04.048]

32. Chen, P.; Xiong, Z.; Xu, J.; Liu, M. Simulation and parameter optimization of high moisture rice drying on combine harvester before threshing. Comput. Electron. Agric.; 2023; 215, 108451. [DOI: https://dx.doi.org/10.1016/j.compag.2023.108451]

33. Liu, W.; Zhou, Y.; Xu, H.; Fu, J.; Zhang, N.; Xie, G.; Zhang, G. Optimization and experiments of the drum longitudinal axial threshing cylinder with rod tooth for rice. Trans. Chin. Soc. Agric. Eng.; 2023; 39, pp. 34-45. [DOI: https://dx.doi.org/10.11975/j.issn.1002-6819.202305076]

34. Xu, C.; Xu, F.; Tang, H.; Wang, J. Determination of characteristics and establishment of discrete element model for whole rice plant. Agronomy; 2023; 13, 2098. [DOI: https://dx.doi.org/10.3390/agronomy13082098]

35. Wang, Q.; Mao, H.; Li, Q. Modelling and simulation of the grain threshing process based on the discrete element method. Comput. Electron. Agric.; 2020; 178, 105790. [DOI: https://dx.doi.org/10.1016/j.compag.2020.105790]

36. Potyondy, D.O.; Cundall, P.A. A bonded-particle model for rock. Int. J. Rock Mech. Min. Sci.; 2004; 41, pp. 1329-1364. [DOI: https://dx.doi.org/10.1016/j.ijrmms.2004.09.011]

37. Coetzee, C.J.; Lombard, S.G. The destemming of grapes: Experiments and discrete element modelling. Biosyst. Eng.; 2013; 114, pp. 232-248. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2012.12.014]

38. Xu, L.; Li, Y.; Ding, L. Contacting mechanics analysis during impact process between rice and threshing component. Trans. CSAE; 2008; 24, pp. 146-149. [DOI: https://dx.doi.org/10.3901/JME.2008.09.177]

39. Li, H.; Li, Y.; Gao, F.; Zhao, Z.; Xu, L. CFD–DEM simulation of material motion in air-and-screen cleaning device. Comput. Electron. Agric.; 2012; 88, pp. 111-119. [DOI: https://dx.doi.org/10.1016/j.compag.2012.07.006]

40. Horabik, J.; Molenda, M. Parameters and contact models for DEM simulations of agricultural granular materials: A review. Biosyst. Eng.; 2016; 147, pp. 206-225. [DOI: https://dx.doi.org/10.1016/j.biosystemseng.2016.02.017]

41. Spagnuolo, S.; Rinaldi, Z.; Donnini, J.; Nanni, A. Physical, mechanical and durability properties of GFRP bars with modified acrylic resin (modar) matrix. Compos. Struct.; 2021; 262, 113557. [DOI: https://dx.doi.org/10.1016/j.compstruct.2021.113557]

42. Stigh, U.; Biel, A. Effects of strain rate on the cohesive properties and fracture process of a pressure sensitive adhesive. Eng. Fract. Mech.; 2018; 203, pp. 266-275. [DOI: https://dx.doi.org/10.1016/j.engfracmech.2018.07.011]

43. Han, D.; Li, W.; Cui, T.; Wang, Q.; Tang, C.; Chen, L.; Xu, L. Biomechanical features and parametric calibration of a bilayer bonded particle model for the cornstalk during harvest. Sci. Rep.; 2024; 14, 20198. [DOI: https://dx.doi.org/10.1038/s41598-024-71005-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/39215138]

44. Sun, K.; Yu, J.; Zhao, J.; Liang, L.; Wang, Y.; Yu, Y. A DEM-based general modeling method and experimental verification for wheat plants in the mature period. Comput. Electron. Agric.; 2023; 214, 108283. [DOI: https://dx.doi.org/10.1016/j.compag.2023.108283]

45. Ashtiani Araghi, H.; Sadeghi, M.; Hemmat, A. Physical properties of two rough rice varieties affected by moisture content. Int. Agrophys.; 2010; 24, pp. 205-207.

46. Wang, Q.; Mao, H.; Li, Q. Simulation of Vibration Response of Flexible Crop Stem Based on Discrete Element Method. Trans. Chin. Soc. Agric. Mach.; 2020; 51, pp. 131-137. [DOI: https://dx.doi.org/10.6041/j.issn.1000-1298.2020.11.014]